Problem 1

Remember the tomato data set generated by Pepe; we first looked at this when we were working on Chapter 9. 35 accessions for seven species were grown in sun and shade.

Assess whether there is evidence for total length (“totleng”) response to shade and whether this response varies by species. Consider whether including accession (“acs”) or shelf (“shelf”) using adaptive priors improves the model fit.

get the data and format

tomato <- read.csv("Assignment_Chapter_09/TomatoR2CSHL.csv")
names(tomato)
##  [1] "shelf"    "flat"     "col"      "row"      "acs"      "trt"     
##  [7] "days"     "date"     "hyp"      "int1"     "int2"     "int3"    
## [13] "int4"     "intleng"  "totleng"  "petleng"  "leafleng" "leafwid" 
## [19] "leafnum"  "ndvi"     "lat"      "lon"      "alt"      "species" 
## [25] "who"
tomato$species_id <- coerce_index(tomato$species)
tomato$shade <- ifelse(tomato$trt=="L",1,0)
tomato$shelf_id <- coerce_index(tomato$shelf)
tomato$species_id <- coerce_index(tomato$species)
tomato$acs_id <- coerce_index(tomato$acs)
tomato.small <- tomato[,c("shelf_id","acs_id","shade","totleng","species_id")]
summary(tomato.small)
##     shelf_id         acs_id         shade           totleng      
##  Min.   :1.000   Min.   : 1.0   Min.   :0.0000   Min.   : 13.59  
##  1st Qu.:2.000   1st Qu.: 9.0   1st Qu.:0.0000   1st Qu.: 39.25  
##  Median :3.000   Median :18.0   Median :1.0000   Median : 50.98  
##  Mean   :3.512   Mean   :18.1   Mean   :0.5089   Mean   : 53.70  
##  3rd Qu.:5.000   3rd Qu.:27.0   3rd Qu.:1.0000   3rd Qu.: 64.76  
##  Max.   :6.000   Max.   :36.0   Max.   :1.0000   Max.   :129.43  
##    species_id   
##  Min.   :1.000  
##  1st Qu.:2.000  
##  Median :3.000  
##  Mean   :2.927  
##  3rd Qu.:4.000  
##  Max.   :5.000

shade effect

mean(tomato.small$totleng[tomato.small$shade==0])
## [1] 42.0596
m.shade <- map2stan(
  alist(
    totleng ~ dnorm(mu,sigma),
    mu <- alpha + 
      beta_shade*shade,
    alpha ~ dnorm(42,10),
    beta_shade ~ dnorm(0,20),
    sigma ~ dcauchy(0,1)
  ),
  data = tomato.small,
  chains=4
)
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1, Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1, Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1, Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1, Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%]  (Sampling)
##  Elapsed Time: 0.364957 seconds (Warm-up)
##                0.304377 seconds (Sampling)
##                0.669334 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 2).
## 
## Chain 2, Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2, Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2, Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2, Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2, Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2, Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2, Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2, Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2, Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2, Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2, Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2, Iteration: 2000 / 2000 [100%]  (Sampling)
##  Elapsed Time: 0.4212 seconds (Warm-up)
##                0.293224 seconds (Sampling)
##                0.714424 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 3).
## 
## Chain 3, Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3, Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3, Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3, Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3, Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3, Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3, Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3, Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3, Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3, Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3, Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3, Iteration: 2000 / 2000 [100%]  (Sampling)
##  Elapsed Time: 0.300826 seconds (Warm-up)
##                0.244408 seconds (Sampling)
##                0.545234 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 4).
## 
## Chain 4, Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4, Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4, Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4, Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4, Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4, Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4, Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4, Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4, Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4, Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4, Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4, Iteration: 2000 / 2000 [100%]  (Sampling)
##  Elapsed Time: 0.274527 seconds (Warm-up)
##                0.257516 seconds (Sampling)
##                0.532043 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 4e-06 seconds (Warm-up)
##                0.000122 seconds (Sampling)
##                0.000126 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
plot(m.shade)

pairs(m.shade)

precis(m.shade)
##             Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha      42.08   0.69      40.95      43.15  1803    1
## beta_shade 22.81   0.95      21.30      24.33  1862    1
## sigma      16.05   0.34      15.51      16.59  2627    1

shade and species

m.shade.species <- map2stan(
  alist(
    totleng ~ dnorm(mu,sigma),
    mu <- alpha + 
      beta_sp[species_id] + 
      beta_shade*shade,
    alpha ~ dnorm(42,10),
    beta_sp[species_id] ~ dnorm(0,10),
    beta_shade ~ dnorm(0,20),
    sigma ~ dcauchy(0,1)
  ),
  data = tomato.small,
  chains=4
)
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## 
## Chain 1, Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 1, Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 1, Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 1, Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 1, Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 1, Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 1, Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1, Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 1, Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1, Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1, Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1, Iteration: 2000 / 2000 [100%]  (Sampling)
##  Elapsed Time: 2.65907 seconds (Warm-up)
##                1.81048 seconds (Sampling)
##                4.46955 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 2).
## 
## Chain 2, Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 2, Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 2, Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 2, Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 2, Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 2, Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 2, Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 2, Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 2, Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 2, Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2, Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2, Iteration: 2000 / 2000 [100%]  (Sampling)
##  Elapsed Time: 2.69747 seconds (Warm-up)
##                1.57889 seconds (Sampling)
##                4.27636 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 3).
## 
## Chain 3, Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 3, Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 3, Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 3, Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 3, Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 3, Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 3, Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 3, Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3, Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3, Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3, Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3, Iteration: 2000 / 2000 [100%]  (Sampling)
##  Elapsed Time: 2.71582 seconds (Warm-up)
##                1.84126 seconds (Sampling)
##                4.55708 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 4).
## 
## Chain 4, Iteration:    1 / 2000 [  0%]  (Warmup)
## Chain 4, Iteration:  200 / 2000 [ 10%]  (Warmup)
## Chain 4, Iteration:  400 / 2000 [ 20%]  (Warmup)
## Chain 4, Iteration:  600 / 2000 [ 30%]  (Warmup)
## Chain 4, Iteration:  800 / 2000 [ 40%]  (Warmup)
## Chain 4, Iteration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4, Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 4, Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 4, Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4, Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4, Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4, Iteration: 2000 / 2000 [100%]  (Sampling)
##  Elapsed Time: 2.66713 seconds (Warm-up)
##                1.65344 seconds (Sampling)
##                4.32057 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 5e-06 seconds (Warm-up)
##                0.000182 seconds (Sampling)
##                0.000187 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
plot(m.shade.species)
pairs(m.shade.species)

precis(m.shade.species,depth = 2)
##             Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha      41.13   4.21      34.38      47.85   641 1.01
## beta_sp[1] -1.29   4.25      -8.31       5.34   444 1.00
## beta_sp[2]  0.12   4.27      -7.08       6.52   659 1.00
## beta_sp[3]  2.42   4.25      -4.12       9.43   663 1.00
## beta_sp[4] -9.41   4.29     -16.71      -2.86   675 1.00
## beta_sp[5]  8.35   4.24       1.44      15.05   667 1.00
## beta_shade 22.99   0.93      21.56      24.48  1757 1.00
## sigma      15.17   0.33      14.66      15.69  1947 1.00
compare(m.shade,m.shade.species) 
##                   WAIC pWAIC dWAIC weight    SE   dSE
## m.shade.species 8348.1   6.9   0.0      1 53.11    NA
## m.shade         8458.4   3.2 110.3      0 52.47 21.97
levels(tomato$species)
## [1] "S. chilense"     "S. chmielewskii" "S. habrochaites" "S. pennellii"   
## [5] "S. peruvianum"

Good evidence for species being important

consider species X shade interaction

m.shade.species.int <- map2stan(
  alist(
    totleng ~ dnorm(mu,sigma),
    mu <- alpha + 
      beta_sp[species_id] + 
      beta_shade*shade + 
      beta_sp_shade[species_id]*shade,
    alpha ~ dnorm(53,20),
    beta_sp[species_id] ~ dnorm(0,10),
    beta_shade ~ dnorm(0,20),
    beta_sp_shade[species_id] ~ dnorm(0,10),
    sigma ~ dcauchy(0,1)
  ),
  data = tomato.small,
  chains=4,
  cores=2
)
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 8e-06 seconds (Warm-up)
##                0.000484 seconds (Sampling)
##                0.000492 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
plot(m.shade.species.int)
par(mfrow=c(1,1))

precis(m.shade.species.int,depth=2)
##                    Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha             41.58   4.42      34.51      48.32   990    1
## beta_sp[1]        -4.24   4.56     -11.35       2.86  1049    1
## beta_sp[2]         2.67   4.54      -3.99      10.04  1044    1
## beta_sp[3]         3.15   4.56      -4.02      10.11  1040    1
## beta_sp[4]       -11.83   4.61     -19.03      -4.53  1058    1
## beta_sp[5]         7.03   4.55      -0.15      14.12  1081    1
## beta_shade        22.19   4.27      15.63      29.17  1075    1
## beta_sp_shade[1]   5.84   4.63      -1.23      13.56  1219    1
## beta_sp_shade[2]  -5.10   4.47     -12.07       2.09  1224    1
## beta_sp_shade[3]  -1.71   4.53      -9.20       5.49  1180    1
## beta_sp_shade[4]   4.43   4.68      -2.83      12.31  1344    1
## beta_sp_shade[5]   2.51   4.57      -4.65       9.80  1194    1
## sigma             15.03   0.33      14.51      15.53  2566    1
levels(tomato$species)
## [1] "S. chilense"     "S. chmielewskii" "S. habrochaites" "S. pennellii"   
## [5] "S. peruvianum"
compare(m.shade,m.shade.species,m.shade.species.int)
##                       WAIC pWAIC dWAIC weight    SE   dSE
## m.shade.species.int 8334.6  10.6   0.0      1 53.30    NA
## m.shade.species     8348.1   6.9  13.5      0 53.11  8.16
## m.shade             8458.4   3.2 123.8      0 52.47 22.14
plot(coeftab(m.shade,m.shade.species,m.shade.species.int))

Reasonable evidence for the interaction, although the 95% posterior intervals overlap among all of the species_shade interactions

Add adaptive prior for accessions

m.shade.species.int.acs <- map2stan(
  alist(
    totleng ~ dnorm(mu,sigma),
    mu <- alpha + 
      alpha_acs[acs_id] +
      beta_sp[species_id] + 
      beta_shade*shade + 
      beta_sp_shade[species_id]*shade,
    alpha ~ dnorm(53,20),
    alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
    beta_sp[species_id] ~ dnorm(0,10),
    beta_shade ~ dnorm(0,20),
    beta_sp_shade[species_id] ~ dnorm(0,10),
    c(sigma,sigma_acs) ~ dcauchy(0,1)
  ),
  data = tomato.small,
  iter = 4000, # Rhat a bit high with defaults
  warmup = 1000,
  chains=4,
  cores=2
)
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 3e-06 seconds (Warm-up)
##                0.000359 seconds (Sampling)
##                0.000362 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]
plot(m.shade.species.int.acs,ask=FALSE)
## Waiting to draw page 2 of 4

## Waiting to draw page 3 of 4

## Waiting to draw page 4 of 4

par(mfrow=c(1,1))

precis(m.shade.species.int.acs,depth = 2)
##                    Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha             41.50   4.44      34.56      48.78  2782    1
## alpha_acs[1]      -5.40   3.68     -11.19       0.52  4841    1
## alpha_acs[2]     -10.63   3.53     -16.48      -5.21  4319    1
## alpha_acs[3]      -3.54   3.42      -9.01       1.83  4277    1
## alpha_acs[4]       3.15   3.38      -1.95       8.75  4305    1
## alpha_acs[5]      13.56   3.41       8.13      18.93  4221    1
## alpha_acs[6]       2.48   3.44      -2.97       7.87  4133    1
## alpha_acs[7]       9.86   3.46       4.30      15.31  4518    1
## alpha_acs[8]       6.35   3.33       1.24      11.84  3916    1
## alpha_acs[9]       3.44   3.37      -1.59       9.08  3752    1
## alpha_acs[10]      3.86   3.90      -2.12      10.29  5649    1
## alpha_acs[11]      2.31   3.32      -2.90       7.74  3848    1
## alpha_acs[12]     -9.78   3.50     -15.44      -4.35  4320    1
## alpha_acs[13]     -4.01   3.66      -9.69       1.93  4923    1
## alpha_acs[14]     -0.69   3.53      -6.04       5.22  4108    1
## alpha_acs[15]     -3.50   4.99     -11.56       4.29  7840    1
## alpha_acs[16]      1.36   3.29      -3.76       6.70  4028    1
## alpha_acs[17]      4.78   3.39      -0.61      10.21  4564    1
## alpha_acs[18]    -12.07   3.53     -17.48      -6.29  4835    1
## alpha_acs[19]     15.04   3.46       9.71      20.77  4066    1
## alpha_acs[20]     -6.18   3.31     -11.25      -0.70  3712    1
## alpha_acs[21]      0.96   3.36      -4.43       6.34  3990    1
## alpha_acs[22]      6.18   4.52      -1.01      13.33  6243    1
## alpha_acs[23]     -2.24   3.68      -8.12       3.56  4884    1
## alpha_acs[24]     -4.28   3.38      -9.49       1.29  4363    1
## alpha_acs[25]     -5.10   3.31     -10.27       0.32  3994    1
## alpha_acs[26]     -0.31   3.55      -5.93       5.42  4978    1
## alpha_acs[27]      1.83   3.33      -3.51       7.01  4483    1
## alpha_acs[28]     -4.60   3.47     -10.44       0.68  4768    1
## alpha_acs[29]      1.14   3.55      -4.56       6.74  4772    1
## alpha_acs[30]      6.78   3.42       1.31      12.18  4699    1
## alpha_acs[31]     10.34   3.49       4.88      15.93  3991    1
## alpha_acs[32]      0.51   3.59      -5.26       6.08  4599    1
## alpha_acs[33]     -3.30   3.54      -8.88       2.39  4628    1
## alpha_acs[34]     -8.61   3.40     -14.11      -3.28  4360    1
## alpha_acs[35]     -5.57   3.53     -11.14      -0.03  4579    1
## alpha_acs[36]     -7.25   4.36     -14.21      -0.35  5457    1
## beta_sp[1]        -4.01   5.00     -12.29       3.59  2889    1
## beta_sp[2]         2.54   4.98      -5.24      10.57  2864    1
## beta_sp[3]         2.65   4.93      -5.24      10.52  2974    1
## beta_sp[4]       -10.69   5.09     -18.66      -2.42  3197    1
## beta_sp[5]         6.23   5.04      -2.03      14.05  2759    1
## beta_shade        22.05   4.43      14.73      28.86  2419    1
## beta_sp_shade[1]   5.72   4.66      -2.02      12.86  2695    1
## beta_sp_shade[2]  -5.12   4.64     -11.88       2.93  2648    1
## beta_sp_shade[3]  -1.89   4.65      -8.95       5.86  2676    1
## beta_sp_shade[4]   5.11   4.80      -2.46      12.78  2827    1
## beta_sp_shade[5]   2.39   4.65      -4.59      10.31  2744    1
## sigma             13.34   0.31      12.85      13.82  8622    1
## sigma_acs          7.47   1.08       5.70       9.05  5679    1
compare(m.shade.species.int,m.shade.species.int.acs)
##                           WAIC pWAIC dWAIC weight    SE   dSE
## m.shade.species.int.acs 8124.3  38.7   0.0      1 55.22    NA
## m.shade.species.int     8334.6  10.6 210.3      0 53.30 28.15
coeftab(m.shade.species.int,m.shade.species.int.acs)
##                  m.shade.species.int m.shade.species.int.acs
## alpha              41.58               41.50                
## beta_sp[1]         -4.24               -4.01                
## beta_sp[2]          2.67                2.54                
## beta_sp[3]          3.15                2.65                
## beta_sp[4]        -11.83              -10.69                
## beta_sp[5]          7.03                6.23                
## beta_shade         22.19               22.05                
## beta_sp_shade[1]    5.84                5.72                
## beta_sp_shade[2]   -5.10               -5.12                
## beta_sp_shade[3]   -1.71               -1.89                
## beta_sp_shade[4]    4.43                5.11                
## beta_sp_shade[5]    2.51                2.39                
## sigma              15.03               13.34                
## alpha_acs[1]          NA                -5.4                
## alpha_acs[2]          NA              -10.63                
## alpha_acs[3]          NA               -3.54                
## alpha_acs[4]          NA                3.15                
## alpha_acs[5]          NA               13.56                
## alpha_acs[6]          NA                2.48                
## alpha_acs[7]          NA                9.86                
## alpha_acs[8]          NA                6.35                
## alpha_acs[9]          NA                3.44                
## alpha_acs[10]         NA                3.86                
## alpha_acs[11]         NA                2.31                
## alpha_acs[12]         NA               -9.78                
## alpha_acs[13]         NA               -4.01                
## alpha_acs[14]         NA               -0.69                
## alpha_acs[15]         NA                -3.5                
## alpha_acs[16]         NA                1.36                
## alpha_acs[17]         NA                4.78                
## alpha_acs[18]         NA              -12.07                
## alpha_acs[19]         NA               15.04                
## alpha_acs[20]         NA               -6.18                
## alpha_acs[21]         NA                0.96                
## alpha_acs[22]         NA                6.18                
## alpha_acs[23]         NA               -2.24                
## alpha_acs[24]         NA               -4.28                
## alpha_acs[25]         NA                -5.1                
## alpha_acs[26]         NA               -0.31                
## alpha_acs[27]         NA                1.83                
## alpha_acs[28]         NA                -4.6                
## alpha_acs[29]         NA                1.14                
## alpha_acs[30]         NA                6.78                
## alpha_acs[31]         NA               10.34                
## alpha_acs[32]         NA                0.51                
## alpha_acs[33]         NA                -3.3                
## alpha_acs[34]         NA               -8.61                
## alpha_acs[35]         NA               -5.57                
## alpha_acs[36]         NA               -7.25                
## sigma_acs             NA                7.47                
## nobs                1008                1008
plot(coeftab(m.shade.species.int,m.shade.species.int.acs),pars=1:13)

Adding acs does not change most of the shared coefficients by much, although it drops sigma.

Add shelf random effect

table(tomato.small$shade, tomato.small$shelf_id)
##    
##       1   2   3   4   5   6
##   0   0   0   0 174 125 196
##   1 161 174 178   0   0   0
table(tomato.small$acs_id,tomato.small$shelf_id)
##     
##       1  2  3  4  5  6
##   1   0  5  5  0  4  6
##   2   2  6  6  3  2  6
##   3   6  6  6  6  4  6
##   4   6  6  6  6  3  6
##   5   6  6  6  5  5  6
##   6   5  5  5  5  5  6
##   7   5  6  6  6  2  6
##   8   7  4  6  9  5  6
##   9   6  5  6  6  5  6
##   10  0  6  2  0  3  5
##   11  6  6  6  9  4  5
##   12  2  4  6  3  5  6
##   13  3  6  5  1  3  6
##   14  2  2  6  3  4  6
##   15  0  2  0  0  0  3
##   16  8  6  6  9  5  6
##   17  6  5  5  7  4  6
##   18  4  6  3  5  4  6
##   19  5  6  5  5  3  6
##   20 10  2  6 11  3  5
##   21  8  4  4  9  5  6
##   22  0  5  0  0  0  3
##   23  3  4  5  0  3  6
##   24  6  6  6  6  5  6
##   25  8  6  6  9  4  6
##   26  3  3  6  4  2  4
##   27  7  6  5  9  4  6
##   28  4  6  5  5  4  5
##   29  4  4  5  5  3  5
##   30  5  6  7  5  3  6
##   31  6  3  6  4  4  5
##   32  3  4  4  6  3  5
##   33  5  5  6  5  1  6
##   34  5  6  5  5  4  6
##   35  3  5  3  3  4  6
##   36  2  1  3  0  3  1

most accessions are on most shelves

m.shade.species.int.acs.shelf <- map2stan(
  alist(
    totleng ~ dnorm(mu,sigma),
    mu <- alpha + 
      alpha_acs[acs_id] +
      alpha_shelf[shelf_id]  +
      beta_sp[species_id] + 
      beta_shade*shade + 
      beta_sp_shade[species_id]*shade,
    alpha ~ dnorm(53,20),
    alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
    alpha_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
    beta_sp[species_id] ~ dnorm(0,10),
    beta_shade ~ dnorm(0,20),
    beta_sp_shade[species_id] ~ dnorm(0,10),
    c(sigma,sigma_acs) ~ dcauchy(0,1),
    sigma_shelf ~ dexp(1) # with only 3 shelves per trt might be hard to estimate w. cauchy
  ),
  data = tomato.small,
  iter = 4000, # Rhat a bit high with defaults
  warmup = 1000,
  chains=4,
  cores=2
)
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 7e-06 seconds (Warm-up)
##                0.000508 seconds (Sampling)
##                0.000515 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]
save.image(file="tomato1212.Rdata")
plot(m.shade.species.int.acs.shelf,ask=FALSE)
## Waiting to draw page 2 of 4

## Waiting to draw page 3 of 4

## Waiting to draw page 4 of 4

par(mfrow=c(1,1))

precis(m.shade.species.int.acs.shelf,depth = 2)
##                    Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha             41.28   4.71      33.63      48.64  3876    1
## alpha_acs[1]      -6.46   3.66     -12.32      -0.60  6404    1
## alpha_acs[2]     -11.35   3.56     -16.99      -5.65  5773    1
## alpha_acs[3]      -3.32   3.37      -8.74       1.94  5714    1
## alpha_acs[4]       3.43   3.41      -2.09       8.70  6028    1
## alpha_acs[5]      13.88   3.42       8.44      19.30  5984    1
## alpha_acs[6]       2.73   3.41      -2.48       8.33  6098    1
## alpha_acs[7]       9.93   3.44       4.42      15.36  5701    1
## alpha_acs[8]       6.89   3.30       1.68      12.23  4995    1
## alpha_acs[9]       3.63   3.34      -1.91       8.76  5073    1
## alpha_acs[10]      2.97   3.91      -2.85       9.52  7398    1
## alpha_acs[11]      2.70   3.34      -2.67       7.97  5772    1
## alpha_acs[12]    -10.18   3.47     -15.85      -4.70  5454    1
## alpha_acs[13]     -4.09   3.66      -9.86       1.74  6426    1
## alpha_acs[14]     -0.79   3.55      -6.38       4.90  5565    1
## alpha_acs[15]     -5.53   4.92     -12.99       2.72  7819    1
## alpha_acs[16]      1.86   3.33      -3.30       7.17  5622    1
## alpha_acs[17]      5.00   3.43      -0.51      10.43  5559    1
## alpha_acs[18]    -12.38   3.55     -17.90      -6.55  5381    1
## alpha_acs[19]     15.11   3.48       9.61      20.59  5765    1
## alpha_acs[20]     -5.13   3.32     -10.55       0.06  5262    1
## alpha_acs[21]      1.60   3.30      -3.79       6.74  5078    1
## alpha_acs[22]      4.56   4.52      -2.25      12.24  9067    1
## alpha_acs[23]     -2.21   3.74      -8.03       3.90  6604    1
## alpha_acs[24]     -4.02   3.35      -9.33       1.34  5907    1
## alpha_acs[25]     -4.56   3.30      -9.54       0.92  5635    1
## alpha_acs[26]     -0.19   3.63      -5.65       6.00  5843    1
## alpha_acs[27]      2.19   3.38      -2.99       7.75  5346    1
## alpha_acs[28]     -4.73   3.48     -10.31       0.86  5552    1
## alpha_acs[29]      1.26   3.55      -4.22       7.07  5770    1
## alpha_acs[30]      6.76   3.44       1.18      12.09  5769    1
## alpha_acs[31]     10.76   3.48       5.45      16.47  5634    1
## alpha_acs[32]      1.30   3.70      -4.34       7.50  6624    1
## alpha_acs[33]     -2.68   3.58      -8.22       3.03  6622    1
## alpha_acs[34]     -8.70   3.45     -14.13      -3.17  6123    1
## alpha_acs[35]     -6.03   3.56     -11.56      -0.18  6366    1
## alpha_acs[36]     -6.46   4.32     -13.46       0.27  7156    1
## alpha_shelf[1]    -3.03   1.88      -6.09      -0.25  5866    1
## alpha_shelf[2]     3.68   1.92       0.71       6.60  5938    1
## alpha_shelf[3]    -0.09   1.87      -2.93       2.86  5612    1
## alpha_shelf[4]    -2.40   1.92      -5.46       0.57  5251    1
## alpha_shelf[5]    -1.09   1.94      -3.96       2.09  5320    1
## alpha_shelf[6]     2.70   1.89      -0.14       5.81  5229    1
## beta_sp[1]        -3.68   4.99     -11.73       4.19  3941    1
## beta_sp[2]         2.65   4.98      -5.19      10.63  3902    1
## beta_sp[3]         2.80   4.96      -4.82      10.94  3488    1
## beta_sp[4]       -11.30   5.03     -19.47      -3.37  4248    1
## beta_sp[5]         6.52   4.99      -1.29      14.46  4143    1
## beta_shade        22.02   5.01      14.11      30.14  3920    1
## beta_sp_shade[1]   5.30   4.61      -2.19      12.54  4287    1
## beta_sp_shade[2]  -5.43   4.57     -12.36       2.22  4070    1
## beta_sp_shade[3]  -1.81   4.56      -9.18       5.31  4081    1
## beta_sp_shade[4]   5.00   4.72      -2.68      12.50  4382    1
## beta_sp_shade[5]   1.92   4.59      -5.84       8.89  4272    1
## sigma             13.05   0.30      12.59      13.53  9577    1
## sigma_acs          7.57   1.10       5.83       9.22  6869    1
## sigma_shelf        2.86   0.91       1.53       4.18  5378    1
compare(m.shade,m.shade.species,m.shade.species.int,m.shade.species.int.acs,m.shade.species.int.acs.shelf)
##                                 WAIC pWAIC dWAIC weight    SE   dSE
## m.shade.species.int.acs.shelf 8083.3  42.3   0.0      1 56.39    NA
## m.shade.species.int.acs       8124.3  38.7  41.0      0 55.22 12.04
## m.shade.species.int           8334.6  10.6 251.3      0 53.30 30.23
## m.shade.species               8348.1   6.9 264.8      0 53.11 31.20
## m.shade                       8458.4   3.2 375.1      0 52.47 38.45
coeftab(m.shade.species.int.acs,m.shade.species.int.acs.shelf)
##                  m.shade.species.int.acs m.shade.species.int.acs.shelf
## alpha              41.50                   41.28                      
## alpha_acs[1]       -5.40                   -6.46                      
## alpha_acs[2]      -10.63                  -11.35                      
## alpha_acs[3]       -3.54                   -3.32                      
## alpha_acs[4]        3.15                    3.43                      
## alpha_acs[5]       13.56                   13.88                      
## alpha_acs[6]        2.48                    2.73                      
## alpha_acs[7]        9.86                    9.93                      
## alpha_acs[8]        6.35                    6.89                      
## alpha_acs[9]        3.44                    3.63                      
## alpha_acs[10]       3.86                    2.97                      
## alpha_acs[11]       2.31                    2.70                      
## alpha_acs[12]      -9.78                  -10.18                      
## alpha_acs[13]      -4.01                   -4.09                      
## alpha_acs[14]      -0.69                   -0.79                      
## alpha_acs[15]      -3.50                   -5.53                      
## alpha_acs[16]       1.36                    1.86                      
## alpha_acs[17]       4.78                    5.00                      
## alpha_acs[18]     -12.07                  -12.38                      
## alpha_acs[19]      15.04                   15.11                      
## alpha_acs[20]      -6.18                   -5.13                      
## alpha_acs[21]       0.96                    1.60                      
## alpha_acs[22]       6.18                    4.56                      
## alpha_acs[23]      -2.24                   -2.21                      
## alpha_acs[24]      -4.28                   -4.02                      
## alpha_acs[25]      -5.10                   -4.56                      
## alpha_acs[26]      -0.31                   -0.19                      
## alpha_acs[27]       1.83                    2.19                      
## alpha_acs[28]      -4.60                   -4.73                      
## alpha_acs[29]       1.14                    1.26                      
## alpha_acs[30]       6.78                    6.76                      
## alpha_acs[31]      10.34                   10.76                      
## alpha_acs[32]       0.51                    1.30                      
## alpha_acs[33]      -3.30                   -2.68                      
## alpha_acs[34]      -8.61                   -8.70                      
## alpha_acs[35]      -5.57                   -6.03                      
## alpha_acs[36]      -7.25                   -6.46                      
## beta_sp[1]         -4.01                   -3.68                      
## beta_sp[2]          2.54                    2.65                      
## beta_sp[3]          2.65                    2.80                      
## beta_sp[4]        -10.69                  -11.30                      
## beta_sp[5]          6.23                    6.52                      
## beta_shade         22.05                   22.02                      
## beta_sp_shade[1]    5.72                    5.30                      
## beta_sp_shade[2]   -5.12                   -5.43                      
## beta_sp_shade[3]   -1.89                   -1.81                      
## beta_sp_shade[4]    5.11                    5.00                      
## beta_sp_shade[5]    2.39                    1.92                      
## sigma              13.34                   13.05                      
## sigma_acs           7.47                    7.57                      
## alpha_shelf[1]        NA                   -3.03                      
## alpha_shelf[2]        NA                    3.68                      
## alpha_shelf[3]        NA                   -0.09                      
## alpha_shelf[4]        NA                    -2.4                      
## alpha_shelf[5]        NA                   -1.09                      
## alpha_shelf[6]        NA                     2.7                      
## sigma_shelf           NA                    2.86                      
## nobs                1008                    1008
plot(coeftab(m.shade.species.int.acs,m.shade.species.int.acs.shelf),pars=c(1,38:50))

including shelf does not change the estimates much, although that model is apparently preferred

reparamterize with no “main” shade effect

m.shade.species.sep.acs.shelf <- map2stan(
  alist(
    totleng ~ dnorm(mu,sigma),
    mu <- alpha + 
      alpha_acs[acs_id] +
      alpha_shelf[shelf_id]  +
      beta_sp[species_id] + 
      beta_sp_shade[species_id]*shade,
    alpha ~ dnorm(53,20),
    alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
    alpha_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
    beta_sp[species_id] ~ dnorm(0,10),
    beta_sp_shade[species_id] ~ dnorm(0,20),
    c(sigma,sigma_acs) ~ dcauchy(0,1),
    sigma_shelf ~ dexp(1) # with only 3 shelves per trt might be hard to estimate w. cauchy
  ),
  data = tomato.small,
  iter = 4000, # Rhat a bit high with defaults
  warmup = 1000,
  chains=4,
  cores=4
)
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 4e-06 seconds (Warm-up)
##                0.000306 seconds (Sampling)
##                0.00031 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]
save.image(file="tomato1212.Rdata")
plot(m.shade.species.sep.acs.shelf,ask=FALSE)
## Waiting to draw page 2 of 4

## Waiting to draw page 3 of 4

## Waiting to draw page 4 of 4

par(mfrow=c(1,1))

precis(m.shade.species.sep.acs.shelf,depth = 2)
##                    Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha             42.01   4.90      34.24      49.73  2020    1
## alpha_acs[1]      -6.49   3.76     -12.47      -0.57  3002    1
## alpha_acs[2]     -11.40   3.57     -17.33      -5.94  3545    1
## alpha_acs[3]      -3.35   3.43      -8.67       2.17  2703    1
## alpha_acs[4]       3.41   3.43      -1.91       9.13  2713    1
## alpha_acs[5]      13.88   3.49       8.29      19.32  2924    1
## alpha_acs[6]       2.70   3.50      -2.95       8.11  2860    1
## alpha_acs[7]       9.92   3.47       4.12      15.17  3288    1
## alpha_acs[8]       6.91   3.21       1.92      12.16  3134    1
## alpha_acs[9]       3.66   3.27      -1.75       8.70  3300    1
## alpha_acs[10]      2.87   3.94      -3.23       9.20  3897    1
## alpha_acs[11]      2.66   3.36      -2.51       8.11  3133    1
## alpha_acs[12]    -10.18   3.42     -15.86      -4.98  3582    1
## alpha_acs[13]     -4.14   3.71     -10.25       1.53  3588    1
## alpha_acs[14]     -0.82   3.50      -6.40       4.71  3933    1
## alpha_acs[15]     -5.58   4.95     -13.28       2.52  5734    1
## alpha_acs[16]      1.78   3.35      -3.72       6.97  2886    1
## alpha_acs[17]      4.99   3.46      -0.34      10.70  3087    1
## alpha_acs[18]    -12.46   3.52     -18.18      -7.01  3141    1
## alpha_acs[19]     15.04   3.51       9.30      20.49  2977    1
## alpha_acs[20]     -5.13   3.20     -10.40      -0.11  3344    1
## alpha_acs[21]      1.59   3.23      -3.82       6.57  3227    1
## alpha_acs[22]      4.55   4.56      -2.97      11.44  4865    1
## alpha_acs[23]     -2.26   3.78      -8.21       3.74  3793    1
## alpha_acs[24]     -4.07   3.43      -9.51       1.44  2774    1
## alpha_acs[25]     -4.59   3.35      -9.99       0.61  2713    1
## alpha_acs[26]     -0.26   3.66      -6.21       5.49  3503    1
## alpha_acs[27]      2.12   3.42      -3.42       7.50  2800    1
## alpha_acs[28]     -4.79   3.48     -10.30       0.80  3069    1
## alpha_acs[29]      1.21   3.60      -4.34       7.14  3121    1
## alpha_acs[30]      6.71   3.48       1.15      12.20  3136    1
## alpha_acs[31]     10.75   3.39       5.38      16.27  3440    1
## alpha_acs[32]      1.22   3.69      -4.89       6.88  3554    1
## alpha_acs[33]     -2.74   3.64      -8.46       3.13  3490    1
## alpha_acs[34]     -8.73   3.48     -14.30      -3.24  3253    1
## alpha_acs[35]     -6.04   3.59     -11.70      -0.26  3440    1
## alpha_acs[36]     -6.42   4.39     -13.11       0.83  4647    1
## alpha_shelf[1]    -2.43   1.96      -5.58       0.56  2130    1
## alpha_shelf[2]     4.35   2.06       1.02       7.47  1925    1
## alpha_shelf[3]     0.55   1.98      -2.62       3.61  1936    1
## alpha_shelf[4]    -3.16   2.03      -6.27       0.03  1923    1
## alpha_shelf[5]    -1.82   2.02      -4.86       1.41  2074    1
## alpha_shelf[6]     1.97   1.94      -0.85       5.16  2052    1
## beta_sp[1]        -3.58   5.01     -11.59       4.46  2411    1
## beta_sp[2]         2.84   5.06      -5.14      10.89  2158    1
## beta_sp[3]         2.90   5.03      -4.88      11.12  2198    1
## beta_sp[4]       -11.16   5.20     -19.30      -2.76  2259    1
## beta_sp[5]         6.62   5.14      -1.37      14.84  2173    1
## beta_sp_shade[1]  25.88   3.17      20.71      30.66  1743    1
## beta_sp_shade[2]  14.93   3.11      10.08      19.88  1671    1
## beta_sp_shade[3]  18.64   3.12      13.91      23.65  1691    1
## beta_sp_shade[4]  25.51   3.46      20.00      30.99  2004    1
## beta_sp_shade[5]  22.47   3.14      17.34      27.20  1672    1
## sigma             13.05   0.29      12.55      13.48  8817    1
## sigma_acs          7.57   1.08       5.92       9.25  5285    1
## sigma_shelf        2.99   0.98       1.56       4.38  2973    1
compare(m.shade,m.shade.species,m.shade.species.int,m.shade.species.int.acs,m.shade.species.int.acs.shelf,m.shade.species.sep.acs.shelf)
##                                 WAIC pWAIC dWAIC weight    SE   dSE
## m.shade.species.int.acs.shelf 8083.3  42.3   0.0   0.55 56.39    NA
## m.shade.species.sep.acs.shelf 8083.7  42.5   0.4   0.45 56.40  0.42
## m.shade.species.int.acs       8124.3  38.7  41.0   0.00 55.22 12.04
## m.shade.species.int           8334.6  10.6 251.3   0.00 53.30 30.23
## m.shade.species               8348.1   6.9 264.8   0.00 53.11 31.20
## m.shade                       8458.4   3.2 375.1   0.00 52.47 38.45
coeftab(m.shade.species.sep.acs.shelf,m.shade.species.int.acs.shelf)
##                  m.shade.species.sep.acs.shelf
## alpha              42.01                      
## alpha_acs[1]       -6.49                      
## alpha_acs[2]      -11.40                      
## alpha_acs[3]       -3.35                      
## alpha_acs[4]        3.41                      
## alpha_acs[5]       13.88                      
## alpha_acs[6]        2.70                      
## alpha_acs[7]        9.92                      
## alpha_acs[8]        6.91                      
## alpha_acs[9]        3.66                      
## alpha_acs[10]       2.87                      
## alpha_acs[11]       2.66                      
## alpha_acs[12]     -10.18                      
## alpha_acs[13]      -4.14                      
## alpha_acs[14]      -0.82                      
## alpha_acs[15]      -5.58                      
## alpha_acs[16]       1.78                      
## alpha_acs[17]       4.99                      
## alpha_acs[18]     -12.46                      
## alpha_acs[19]      15.04                      
## alpha_acs[20]      -5.13                      
## alpha_acs[21]       1.59                      
## alpha_acs[22]       4.55                      
## alpha_acs[23]      -2.26                      
## alpha_acs[24]      -4.07                      
## alpha_acs[25]      -4.59                      
## alpha_acs[26]      -0.26                      
## alpha_acs[27]       2.12                      
## alpha_acs[28]      -4.79                      
## alpha_acs[29]       1.21                      
## alpha_acs[30]       6.71                      
## alpha_acs[31]      10.75                      
## alpha_acs[32]       1.22                      
## alpha_acs[33]      -2.74                      
## alpha_acs[34]      -8.73                      
## alpha_acs[35]      -6.04                      
## alpha_acs[36]      -6.42                      
## alpha_shelf[1]     -2.43                      
## alpha_shelf[2]      4.35                      
## alpha_shelf[3]      0.55                      
## alpha_shelf[4]     -3.16                      
## alpha_shelf[5]     -1.82                      
## alpha_shelf[6]      1.97                      
## beta_sp[1]         -3.58                      
## beta_sp[2]          2.84                      
## beta_sp[3]           2.9                      
## beta_sp[4]        -11.16                      
## beta_sp[5]          6.62                      
## beta_sp_shade[1]   25.88                      
## beta_sp_shade[2]   14.93                      
## beta_sp_shade[3]   18.64                      
## beta_sp_shade[4]   25.51                      
## beta_sp_shade[5]   22.47                      
## sigma              13.05                      
## sigma_acs           7.57                      
## sigma_shelf         2.99                      
## beta_shade            NA                      
## nobs                1008                      
##                  m.shade.species.int.acs.shelf
## alpha              41.28                      
## alpha_acs[1]       -6.46                      
## alpha_acs[2]      -11.35                      
## alpha_acs[3]       -3.32                      
## alpha_acs[4]        3.43                      
## alpha_acs[5]       13.88                      
## alpha_acs[6]        2.73                      
## alpha_acs[7]        9.93                      
## alpha_acs[8]        6.89                      
## alpha_acs[9]        3.63                      
## alpha_acs[10]       2.97                      
## alpha_acs[11]       2.70                      
## alpha_acs[12]     -10.18                      
## alpha_acs[13]      -4.09                      
## alpha_acs[14]      -0.79                      
## alpha_acs[15]      -5.53                      
## alpha_acs[16]       1.86                      
## alpha_acs[17]       5.00                      
## alpha_acs[18]     -12.38                      
## alpha_acs[19]      15.11                      
## alpha_acs[20]      -5.13                      
## alpha_acs[21]       1.60                      
## alpha_acs[22]       4.56                      
## alpha_acs[23]      -2.21                      
## alpha_acs[24]      -4.02                      
## alpha_acs[25]      -4.56                      
## alpha_acs[26]      -0.19                      
## alpha_acs[27]       2.19                      
## alpha_acs[28]      -4.73                      
## alpha_acs[29]       1.26                      
## alpha_acs[30]       6.76                      
## alpha_acs[31]      10.76                      
## alpha_acs[32]       1.30                      
## alpha_acs[33]      -2.68                      
## alpha_acs[34]      -8.70                      
## alpha_acs[35]      -6.03                      
## alpha_acs[36]      -6.46                      
## alpha_shelf[1]     -3.03                      
## alpha_shelf[2]      3.68                      
## alpha_shelf[3]     -0.09                      
## alpha_shelf[4]     -2.40                      
## alpha_shelf[5]     -1.09                      
## alpha_shelf[6]      2.70                      
## beta_sp[1]         -3.68                      
## beta_sp[2]          2.65                      
## beta_sp[3]           2.8                      
## beta_sp[4]        -11.30                      
## beta_sp[5]          6.52                      
## beta_sp_shade[1]    5.30                      
## beta_sp_shade[2]   -5.43                      
## beta_sp_shade[3]   -1.81                      
## beta_sp_shade[4]    5.00                      
## beta_sp_shade[5]    1.92                      
## sigma              13.05                      
## sigma_acs           7.57                      
## sigma_shelf         2.86                      
## beta_shade         22.02                      
## nobs                1008
plot(coeftab(m.shade.species.sep.acs.shelf))

This parameterization leads to the same WAIC, so it is an equivalent model. Interestingly though, the SD on the shade coefficients is smaller, presumably because one fewer parameter is being estimated.

Can we do no main alpha?

m.shade.species.sep.acs.shelf.no.alpha <- map2stan(
  alist(
    totleng ~ dnorm(mu,sigma),
    mu <- 
      alpha_acs[acs_id] +
      alpha_shelf[shelf_id]  +
      beta_sp[species_id] + 
      beta_sp_shade[species_id]*shade,
    alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
    alpha_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
    beta_sp[species_id] ~ dnorm(42,20),
    beta_sp_shade[species_id] ~ dnorm(0,20),
    c(sigma,sigma_acs) ~ dcauchy(0,1),
    sigma_shelf ~ dexp(1) # with only 3 shelves per trt might be hard to estimate w. cauchy
  ),
  data = tomato.small,
  iter = 4000, # Rhat a bit high with defaults
  warmup = 1000,
  chains=4,
  cores=4
)
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 6e-06 seconds (Warm-up)
##                0.000312 seconds (Sampling)
##                0.000318 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]
save.image(file="tomato1212.Rdata")
plot(m.shade.species.sep.acs.shelf.no.alpha,ask=FALSE)
## Waiting to draw page 2 of 4

## Waiting to draw page 3 of 4

## Waiting to draw page 4 of 4

par(mfrow=c(1,1))

precis(m.shade.species.sep.acs.shelf.no.alpha,depth = 2)
##                    Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha_acs[1]      -6.79   3.75     -12.75      -0.89  3134    1
## alpha_acs[2]     -11.67   3.64     -17.38      -5.74  3279    1
## alpha_acs[3]      -3.60   3.48      -9.29       1.83  2720    1
## alpha_acs[4]       3.13   3.47      -2.59       8.43  2819    1
## alpha_acs[5]      13.62   3.47       7.92      18.95  2685    1
## alpha_acs[6]       2.41   3.48      -3.06       8.08  2735    1
## alpha_acs[7]       9.59   3.52       4.14      15.31  2834    1
## alpha_acs[8]       6.75   3.34       1.18      11.88  3047    1
## alpha_acs[9]       3.53   3.38      -1.99       8.78  3022    1
## alpha_acs[10]      3.75   3.99      -2.58      10.09  3719    1
## alpha_acs[11]      2.34   3.46      -3.11       7.84  2778    1
## alpha_acs[12]    -10.30   3.52     -15.74      -4.51  3340    1
## alpha_acs[13]     -3.29   3.76      -9.11       2.79  3303    1
## alpha_acs[14]     -0.95   3.63      -6.66       4.94  3578    1
## alpha_acs[15]     -5.61   4.93     -13.48       2.27  6084    1
## alpha_acs[16]      1.48   3.39      -3.70       7.11  2573    1
## alpha_acs[17]      5.14   3.48      -0.59      10.54  2873    1
## alpha_acs[18]    -12.26   3.62     -17.81      -6.37  3211    1
## alpha_acs[19]     14.69   3.57       9.26      20.53  2797    1
## alpha_acs[20]     -5.27   3.34     -10.70      -0.03  2946    1
## alpha_acs[21]      1.46   3.38      -3.93       6.84  3113    1
## alpha_acs[22]      5.21   4.57      -2.08      12.54  4562    1
## alpha_acs[23]     -1.44   3.82      -7.39       4.78  3462    1
## alpha_acs[24]     -4.32   3.45      -9.58       1.43  2600    1
## alpha_acs[25]     -4.84   3.42     -10.35       0.59  2649    1
## alpha_acs[26]     -0.09   3.69      -5.94       5.81  3262    1
## alpha_acs[27]      2.28   3.44      -3.45       7.49  2836    1
## alpha_acs[28]     -4.63   3.52     -10.19       1.04  2968    1
## alpha_acs[29]      1.38   3.60      -4.47       6.97  3066    1
## alpha_acs[30]      6.85   3.47       1.28      12.28  2882    1
## alpha_acs[31]     10.66   3.49       5.08      16.15  3252    1
## alpha_acs[32]      2.11   3.71      -3.79       8.02  3079    1
## alpha_acs[33]     -1.90   3.66      -7.62       4.01  3242    1
## alpha_acs[34]     -9.04   3.51     -14.89      -3.72  2867    1
## alpha_acs[35]     -6.39   3.60     -12.04      -0.66  3076    1
## alpha_acs[36]     -5.70   4.38     -12.58       1.37  4809    1
## alpha_shelf[1]    -2.38   1.96      -5.57       0.60  2963    1
## alpha_shelf[2]     4.41   2.06       1.35       7.73  2804    1
## alpha_shelf[3]     0.60   1.97      -2.36       3.78  2811    1
## alpha_shelf[4]    -3.05   1.96      -5.99       0.12  1925    1
## alpha_shelf[5]    -1.68   1.95      -4.77       1.37  2061    1
## alpha_shelf[6]     2.10   1.88      -0.89       5.00  2010    1
## beta_sp[1]        38.09   3.49      32.76      43.85  1965    1
## beta_sp[2]        45.03   3.55      39.32      50.61  1930    1
## beta_sp[3]        44.91   3.38      39.77      50.58  1980    1
## beta_sp[4]        29.62   3.74      23.65      35.58  2207    1
## beta_sp[5]        48.90   3.55      43.22      54.42  1953    1
## beta_sp_shade[1]  26.03   3.08      21.38      31.05  2237    1
## beta_sp_shade[2]  14.96   3.04      10.49      20.03  2169    1
## beta_sp_shade[3]  18.69   3.02      14.13      23.70  2219    1
## beta_sp_shade[4]  25.85   3.34      20.73      31.29  2538    1
## beta_sp_shade[5]  22.44   3.05      17.62      27.21  2178    1
## sigma             13.05   0.30      12.58      13.53 10166    1
## sigma_acs          7.57   1.09       5.81       9.16  5217    1
## sigma_shelf        2.98   0.96       1.57       4.34  3519    1
compare(m.shade,m.shade.species,m.shade.species.int,m.shade.species.int.acs,m.shade.species.int.acs.shelf,m.shade.species.sep.acs.shelf,m.shade.species.sep.acs.shelf.no.alpha)
##                                          WAIC pWAIC dWAIC weight    SE
## m.shade.species.sep.acs.shelf.no.alpha 8083.3  42.3   0.0   0.36 56.42
## m.shade.species.int.acs.shelf          8083.3  42.3   0.0   0.35 56.39
## m.shade.species.sep.acs.shelf          8083.7  42.5   0.5   0.28 56.40
## m.shade.species.int.acs                8124.3  38.7  41.0   0.00 55.22
## m.shade.species.int                    8334.6  10.6 251.4   0.00 53.30
## m.shade.species                        8348.1   6.9 264.9   0.00 53.11
## m.shade                                8458.4   3.2 375.1   0.00 52.47
##                                          dSE
## m.shade.species.sep.acs.shelf.no.alpha    NA
## m.shade.species.int.acs.shelf           0.54
## m.shade.species.sep.acs.shelf           0.32
## m.shade.species.int.acs                12.17
## m.shade.species.int                    30.19
## m.shade.species                        31.19
## m.shade                                38.46
coeftab(m.shade.species.sep.acs.shelf,m.shade.species.int.acs.shelf)
##                  m.shade.species.sep.acs.shelf
## alpha              42.01                      
## alpha_acs[1]       -6.49                      
## alpha_acs[2]      -11.40                      
## alpha_acs[3]       -3.35                      
## alpha_acs[4]        3.41                      
## alpha_acs[5]       13.88                      
## alpha_acs[6]        2.70                      
## alpha_acs[7]        9.92                      
## alpha_acs[8]        6.91                      
## alpha_acs[9]        3.66                      
## alpha_acs[10]       2.87                      
## alpha_acs[11]       2.66                      
## alpha_acs[12]     -10.18                      
## alpha_acs[13]      -4.14                      
## alpha_acs[14]      -0.82                      
## alpha_acs[15]      -5.58                      
## alpha_acs[16]       1.78                      
## alpha_acs[17]       4.99                      
## alpha_acs[18]     -12.46                      
## alpha_acs[19]      15.04                      
## alpha_acs[20]      -5.13                      
## alpha_acs[21]       1.59                      
## alpha_acs[22]       4.55                      
## alpha_acs[23]      -2.26                      
## alpha_acs[24]      -4.07                      
## alpha_acs[25]      -4.59                      
## alpha_acs[26]      -0.26                      
## alpha_acs[27]       2.12                      
## alpha_acs[28]      -4.79                      
## alpha_acs[29]       1.21                      
## alpha_acs[30]       6.71                      
## alpha_acs[31]      10.75                      
## alpha_acs[32]       1.22                      
## alpha_acs[33]      -2.74                      
## alpha_acs[34]      -8.73                      
## alpha_acs[35]      -6.04                      
## alpha_acs[36]      -6.42                      
## alpha_shelf[1]     -2.43                      
## alpha_shelf[2]      4.35                      
## alpha_shelf[3]      0.55                      
## alpha_shelf[4]     -3.16                      
## alpha_shelf[5]     -1.82                      
## alpha_shelf[6]      1.97                      
## beta_sp[1]         -3.58                      
## beta_sp[2]          2.84                      
## beta_sp[3]           2.9                      
## beta_sp[4]        -11.16                      
## beta_sp[5]          6.62                      
## beta_sp_shade[1]   25.88                      
## beta_sp_shade[2]   14.93                      
## beta_sp_shade[3]   18.64                      
## beta_sp_shade[4]   25.51                      
## beta_sp_shade[5]   22.47                      
## sigma              13.05                      
## sigma_acs           7.57                      
## sigma_shelf         2.99                      
## beta_shade            NA                      
## nobs                1008                      
##                  m.shade.species.int.acs.shelf
## alpha              41.28                      
## alpha_acs[1]       -6.46                      
## alpha_acs[2]      -11.35                      
## alpha_acs[3]       -3.32                      
## alpha_acs[4]        3.43                      
## alpha_acs[5]       13.88                      
## alpha_acs[6]        2.73                      
## alpha_acs[7]        9.93                      
## alpha_acs[8]        6.89                      
## alpha_acs[9]        3.63                      
## alpha_acs[10]       2.97                      
## alpha_acs[11]       2.70                      
## alpha_acs[12]     -10.18                      
## alpha_acs[13]      -4.09                      
## alpha_acs[14]      -0.79                      
## alpha_acs[15]      -5.53                      
## alpha_acs[16]       1.86                      
## alpha_acs[17]       5.00                      
## alpha_acs[18]     -12.38                      
## alpha_acs[19]      15.11                      
## alpha_acs[20]      -5.13                      
## alpha_acs[21]       1.60                      
## alpha_acs[22]       4.56                      
## alpha_acs[23]      -2.21                      
## alpha_acs[24]      -4.02                      
## alpha_acs[25]      -4.56                      
## alpha_acs[26]      -0.19                      
## alpha_acs[27]       2.19                      
## alpha_acs[28]      -4.73                      
## alpha_acs[29]       1.26                      
## alpha_acs[30]       6.76                      
## alpha_acs[31]      10.76                      
## alpha_acs[32]       1.30                      
## alpha_acs[33]      -2.68                      
## alpha_acs[34]      -8.70                      
## alpha_acs[35]      -6.03                      
## alpha_acs[36]      -6.46                      
## alpha_shelf[1]     -3.03                      
## alpha_shelf[2]      3.68                      
## alpha_shelf[3]     -0.09                      
## alpha_shelf[4]     -2.40                      
## alpha_shelf[5]     -1.09                      
## alpha_shelf[6]      2.70                      
## beta_sp[1]         -3.68                      
## beta_sp[2]          2.65                      
## beta_sp[3]           2.8                      
## beta_sp[4]        -11.30                      
## beta_sp[5]          6.52                      
## beta_sp_shade[1]    5.30                      
## beta_sp_shade[2]   -5.43                      
## beta_sp_shade[3]   -1.81                      
## beta_sp_shade[4]    5.00                      
## beta_sp_shade[5]    1.92                      
## sigma              13.05                      
## sigma_acs           7.57                      
## sigma_shelf         2.86                      
## beta_shade         22.02                      
## nobs                1008
plot(coeftab(m.shade.species.sep.acs.shelf.no.alpha))

So long as an appropriate prior for the species intercepts is used this is OK. The other reasonable paramaterization would be to have the alpha befor the first species.

Consider different slopes for accessions

m.shade.species.sep.acs.int.shelf.no.alpha <- map2stan(
  alist(
    totleng ~ dnorm(mu,sigma),
    mu <- alpha_acs[acs_id] +
      alpha_shelf[shelf_id]  +
      beta_sp[species_id] + 
      beta_acs_shade[acs_id]*shade,
    alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
    alpha_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
    beta_sp[species_id] ~ dnorm(42,10),
    beta_acs_shade[acs_id] ~ dnorm(20,20),
    c(sigma,sigma_acs,sigma_shelf) ~ dcauchy(0,1)
  ),
  data = tomato.small,
  iter = 4000, # Rhat a bit high with defaults
  warmup = 1000,
  chains=4,
  cores=4
)
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 3e-06 seconds (Warm-up)
##                0.000412 seconds (Sampling)
##                0.000415 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]
save.image("tomato1212.Rdata")
plot(m.shade.species.sep.acs.int.shelf.no.alpha,ask=FALSE)
## Waiting to draw page 2 of 6

## Waiting to draw page 3 of 6

## Waiting to draw page 4 of 6

## Waiting to draw page 5 of 6

## Waiting to draw page 6 of 6

par(mfrow=c(1,1))

precis(m.shade.species.sep.acs.int.shelf.no.alpha,depth = 2)
##                     Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha_acs[1]       -1.88   3.27      -7.08       3.26 12000    1
## alpha_acs[2]       -4.05   3.29      -9.52       1.06 12000    1
## alpha_acs[3]       -4.01   3.02      -8.77       0.88 12000    1
## alpha_acs[4]        2.83   3.05      -1.90       7.77 12000    1
## alpha_acs[5]        8.06   3.17       2.85      12.93 12000    1
## alpha_acs[6]       -0.15   3.01      -5.00       4.57 12000    1
## alpha_acs[7]        5.86   3.18       0.88      10.96 12000    1
## alpha_acs[8]        0.12   2.86      -4.56       4.54 12000    1
## alpha_acs[9]        2.53   2.95      -2.13       7.23 12000    1
## alpha_acs[10]      -0.09   3.42      -5.69       5.19 12000    1
## alpha_acs[11]       3.03   2.96      -1.82       7.65 12000    1
## alpha_acs[12]      -4.98   3.12      -9.93       0.06 12000    1
## alpha_acs[13]       0.25   3.32      -5.26       5.43 12000    1
## alpha_acs[14]       2.08   3.11      -3.05       6.82 12000    1
## alpha_acs[15]      -1.06   3.95      -7.59       4.98 12000    1
## alpha_acs[16]       2.69   2.88      -2.07       7.10 12000    1
## alpha_acs[17]       4.76   3.03       0.22       9.78 12000    1
## alpha_acs[18]      -7.57   3.21     -12.78      -2.63 12000    1
## alpha_acs[19]       5.15   3.16       0.24      10.35 12000    1
## alpha_acs[20]      -3.38   2.89      -8.12       1.01 12000    1
## alpha_acs[21]      -0.07   2.85      -4.51       4.55 12000    1
## alpha_acs[22]       0.55   3.95      -5.38       7.24 12000    1
## alpha_acs[23]      -1.77   3.42      -7.32       3.57 12000    1
## alpha_acs[24]      -2.48   2.96      -7.32       2.11 12000    1
## alpha_acs[25]      -1.83   2.91      -6.49       2.76 12000    1
## alpha_acs[26]      -2.09   3.29      -7.23       3.31 12000    1
## alpha_acs[27]       1.02   2.92      -3.59       5.71 12000    1
## alpha_acs[28]      -0.43   3.12      -5.41       4.51 12000    1
## alpha_acs[29]      -0.57   3.12      -5.44       4.50 12000    1
## alpha_acs[30]       3.83   3.12      -1.21       8.75 12000    1
## alpha_acs[31]       5.51   3.17       0.45      10.55 12000    1
## alpha_acs[32]      -0.14   3.14      -5.00       4.96 12000    1
## alpha_acs[33]       0.40   3.25      -4.93       5.40 12000    1
## alpha_acs[34]      -6.86   3.14     -11.89      -1.88 12000    1
## alpha_acs[35]      -4.02   3.13      -8.91       1.03 12000    1
## alpha_acs[36]      -2.10   3.88      -8.21       4.05 12000    1
## alpha_shelf[1]     -3.10   1.91      -6.14      -0.15  6532    1
## alpha_shelf[2]      4.20   1.92       1.15       7.16  6021    1
## alpha_shelf[3]      0.41   1.88      -2.56       3.33  6135    1
## alpha_shelf[4]     -2.56   1.78      -5.36       0.25  5320    1
## alpha_shelf[5]     -1.56   1.81      -4.18       1.50  5755    1
## alpha_shelf[6]      2.30   1.73      -0.28       5.16  5317    1
## beta_sp[1]         37.82   2.60      33.70      41.95  6643    1
## beta_sp[2]         44.54   2.56      40.49      48.58  7021    1
## beta_sp[3]         45.15   2.53      40.98      48.96  5905    1
## beta_sp[4]         29.80   2.78      25.65      34.53  7823    1
## beta_sp[5]         48.55   2.58      44.44      52.66  6719    1
## beta_acs_shade[1]   5.52   5.34      -3.22      13.76 12000    1
## beta_acs_shade[2]   9.18   4.93       1.15      16.87  8680    1
## beta_acs_shade[3]  17.79   4.39      11.07      25.17 12000    1
## beta_acs_shade[4]  15.83   4.45       8.40      22.67 12000    1
## beta_acs_shade[5]  25.38   4.51      17.87      32.31  8931    1
## beta_acs_shade[6]  21.93   4.58      14.74      29.19 12000    1
## beta_acs_shade[7]  29.12   4.55      22.18      36.64 12000    1
## beta_acs_shade[8]  33.92   4.33      26.84      40.65 12000    1
## beta_acs_shade[9]  19.81   4.40      12.96      26.95 12000    1
## beta_acs_shade[10] 34.62   5.65      25.42      43.49 12000    1
## beta_acs_shade[11] 21.04   4.27      13.82      27.46 12000    1
## beta_acs_shade[12]  7.13   4.89      -0.61      14.91 12000    1
## beta_acs_shade[13] 18.76   4.84      11.32      26.87 12000    1
## beta_acs_shade[14]  9.22   5.10       1.01      17.28 12000    1
## beta_acs_shade[15]  1.91   8.92     -11.61      16.61 12000    1
## beta_acs_shade[16] 20.01   4.14      13.71      26.88  9347    1
## beta_acs_shade[17] 26.18   4.49      19.47      33.80 12000    1
## beta_acs_shade[18] 18.01   4.84      10.34      25.81 12000    1
## beta_acs_shade[19] 41.62   4.64      34.17      48.97 12000    1
## beta_acs_shade[20] 14.98   4.30       8.33      22.08 12000    1
## beta_acs_shade[21] 22.02   4.37      15.24      29.14 12000    1
## beta_acs_shade[22] 35.36   6.95      24.57      46.50 12000    1
## beta_acs_shade[23] 27.11   5.14      19.04      35.33 12000    1
## beta_acs_shade[24] 12.70   4.32       5.95      19.66 12000    1
## beta_acs_shade[25] 10.05   4.16       3.54      16.82  8890    1
## beta_acs_shade[26] 31.63   5.12      23.35      39.70 12000    1
## beta_acs_shade[27] 29.20   4.27      22.37      35.99 12000    1
## beta_acs_shade[28] 17.83   4.64      10.36      25.03 12000    1
## beta_acs_shade[29] 31.18   4.85      23.25      38.76 12000    1
## beta_acs_shade[30] 31.55   4.43      24.31      38.43 12000    1
## beta_acs_shade[31] 27.45   4.63      20.25      35.05 12000    1
## beta_acs_shade[32] 31.17   4.95      23.48      39.18 12000    1
## beta_acs_shade[33] 21.07   4.70      13.45      28.37 12000    1
## beta_acs_shade[34] 20.51   4.64      12.99      27.86 12000    1
## beta_acs_shade[35] 19.18   5.01      11.42      27.42 12000    1
## beta_acs_shade[36] 19.26   6.52       8.91      29.80 12000    1
## sigma              12.73   0.29      12.28      13.21 12000    1
## sigma_acs           4.74   0.91       3.26       6.12  4913    1
## sigma_shelf         3.32   1.32       1.57       4.99  6746    1
compare(m.shade.species.sep.acs.int.shelf.no.alpha,m.shade.species.sep.acs.shelf.no.alpha)
##                                              WAIC pWAIC dWAIC weight    SE
## m.shade.species.sep.acs.int.shelf.no.alpha 8062.2  67.3   0.0      1 58.73
## m.shade.species.sep.acs.shelf.no.alpha     8083.3  42.3  21.1      0 56.42
##                                              dSE
## m.shade.species.sep.acs.int.shelf.no.alpha    NA
## m.shade.species.sep.acs.shelf.no.alpha     20.36
coeftab(m.shade.species.sep.acs.int.shelf.no.alpha,m.shade.species.sep.acs.shelf.no.alpha)
##                    m.shade.species.sep.acs.int.shelf.no.alpha
## alpha_acs[1]         -1.88                                   
## alpha_acs[2]         -4.05                                   
## alpha_acs[3]         -4.01                                   
## alpha_acs[4]          2.83                                   
## alpha_acs[5]          8.06                                   
## alpha_acs[6]         -0.15                                   
## alpha_acs[7]          5.86                                   
## alpha_acs[8]          0.12                                   
## alpha_acs[9]          2.53                                   
## alpha_acs[10]        -0.09                                   
## alpha_acs[11]         3.03                                   
## alpha_acs[12]        -4.98                                   
## alpha_acs[13]         0.25                                   
## alpha_acs[14]         2.08                                   
## alpha_acs[15]        -1.06                                   
## alpha_acs[16]         2.69                                   
## alpha_acs[17]         4.76                                   
## alpha_acs[18]        -7.57                                   
## alpha_acs[19]         5.15                                   
## alpha_acs[20]        -3.38                                   
## alpha_acs[21]        -0.07                                   
## alpha_acs[22]         0.55                                   
## alpha_acs[23]        -1.77                                   
## alpha_acs[24]        -2.48                                   
## alpha_acs[25]        -1.83                                   
## alpha_acs[26]        -2.09                                   
## alpha_acs[27]         1.02                                   
## alpha_acs[28]        -0.43                                   
## alpha_acs[29]        -0.57                                   
## alpha_acs[30]         3.83                                   
## alpha_acs[31]         5.51                                   
## alpha_acs[32]        -0.14                                   
## alpha_acs[33]          0.4                                   
## alpha_acs[34]        -6.86                                   
## alpha_acs[35]        -4.02                                   
## alpha_acs[36]         -2.1                                   
## alpha_shelf[1]       -3.10                                   
## alpha_shelf[2]        4.20                                   
## alpha_shelf[3]        0.41                                   
## alpha_shelf[4]       -2.56                                   
## alpha_shelf[5]       -1.56                                   
## alpha_shelf[6]         2.3                                   
## beta_sp[1]           37.82                                   
## beta_sp[2]           44.54                                   
## beta_sp[3]           45.15                                   
## beta_sp[4]           29.80                                   
## beta_sp[5]           48.55                                   
## beta_acs_shade[1]     5.52                                   
## beta_acs_shade[2]     9.18                                   
## beta_acs_shade[3]    17.79                                   
## beta_acs_shade[4]    15.83                                   
## beta_acs_shade[5]    25.38                                   
## beta_acs_shade[6]    21.93                                   
## beta_acs_shade[7]    29.12                                   
## beta_acs_shade[8]    33.92                                   
## beta_acs_shade[9]    19.81                                   
## beta_acs_shade[10]   34.62                                   
## beta_acs_shade[11]   21.04                                   
## beta_acs_shade[12]    7.13                                   
## beta_acs_shade[13]   18.76                                   
## beta_acs_shade[14]    9.22                                   
## beta_acs_shade[15]    1.91                                   
## beta_acs_shade[16]   20.01                                   
## beta_acs_shade[17]   26.18                                   
## beta_acs_shade[18]   18.01                                   
## beta_acs_shade[19]   41.62                                   
## beta_acs_shade[20]   14.98                                   
## beta_acs_shade[21]   22.02                                   
## beta_acs_shade[22]   35.36                                   
## beta_acs_shade[23]   27.11                                   
## beta_acs_shade[24]    12.7                                   
## beta_acs_shade[25]   10.05                                   
## beta_acs_shade[26]   31.63                                   
## beta_acs_shade[27]    29.2                                   
## beta_acs_shade[28]   17.83                                   
## beta_acs_shade[29]   31.18                                   
## beta_acs_shade[30]   31.55                                   
## beta_acs_shade[31]   27.45                                   
## beta_acs_shade[32]   31.17                                   
## beta_acs_shade[33]   21.07                                   
## beta_acs_shade[34]   20.51                                   
## beta_acs_shade[35]   19.18                                   
## beta_acs_shade[36]   19.26                                   
## sigma                12.73                                   
## sigma_acs             4.74                                   
## sigma_shelf           3.32                                   
## beta_sp_shade[1]        NA                                   
## beta_sp_shade[2]        NA                                   
## beta_sp_shade[3]        NA                                   
## beta_sp_shade[4]        NA                                   
## beta_sp_shade[5]        NA                                   
## nobs                  1008                                   
##                    m.shade.species.sep.acs.shelf.no.alpha
## alpha_acs[1]         -6.79                               
## alpha_acs[2]        -11.67                               
## alpha_acs[3]         -3.60                               
## alpha_acs[4]          3.13                               
## alpha_acs[5]         13.62                               
## alpha_acs[6]          2.41                               
## alpha_acs[7]          9.59                               
## alpha_acs[8]          6.75                               
## alpha_acs[9]          3.53                               
## alpha_acs[10]         3.75                               
## alpha_acs[11]         2.34                               
## alpha_acs[12]       -10.30                               
## alpha_acs[13]        -3.29                               
## alpha_acs[14]        -0.95                               
## alpha_acs[15]        -5.61                               
## alpha_acs[16]         1.48                               
## alpha_acs[17]         5.14                               
## alpha_acs[18]       -12.26                               
## alpha_acs[19]        14.69                               
## alpha_acs[20]        -5.27                               
## alpha_acs[21]         1.46                               
## alpha_acs[22]         5.21                               
## alpha_acs[23]        -1.44                               
## alpha_acs[24]        -4.32                               
## alpha_acs[25]        -4.84                               
## alpha_acs[26]        -0.09                               
## alpha_acs[27]         2.28                               
## alpha_acs[28]        -4.63                               
## alpha_acs[29]         1.38                               
## alpha_acs[30]         6.85                               
## alpha_acs[31]        10.66                               
## alpha_acs[32]         2.11                               
## alpha_acs[33]         -1.9                               
## alpha_acs[34]        -9.04                               
## alpha_acs[35]        -6.39                               
## alpha_acs[36]         -5.7                               
## alpha_shelf[1]       -2.38                               
## alpha_shelf[2]        4.41                               
## alpha_shelf[3]        0.60                               
## alpha_shelf[4]       -3.05                               
## alpha_shelf[5]       -1.68                               
## alpha_shelf[6]         2.1                               
## beta_sp[1]           38.09                               
## beta_sp[2]           45.03                               
## beta_sp[3]           44.91                               
## beta_sp[4]           29.62                               
## beta_sp[5]           48.90                               
## beta_acs_shade[1]       NA                               
## beta_acs_shade[2]       NA                               
## beta_acs_shade[3]       NA                               
## beta_acs_shade[4]       NA                               
## beta_acs_shade[5]       NA                               
## beta_acs_shade[6]       NA                               
## beta_acs_shade[7]       NA                               
## beta_acs_shade[8]       NA                               
## beta_acs_shade[9]       NA                               
## beta_acs_shade[10]      NA                               
## beta_acs_shade[11]      NA                               
## beta_acs_shade[12]      NA                               
## beta_acs_shade[13]      NA                               
## beta_acs_shade[14]      NA                               
## beta_acs_shade[15]      NA                               
## beta_acs_shade[16]      NA                               
## beta_acs_shade[17]      NA                               
## beta_acs_shade[18]      NA                               
## beta_acs_shade[19]      NA                               
## beta_acs_shade[20]      NA                               
## beta_acs_shade[21]      NA                               
## beta_acs_shade[22]      NA                               
## beta_acs_shade[23]      NA                               
## beta_acs_shade[24]      NA                               
## beta_acs_shade[25]      NA                               
## beta_acs_shade[26]      NA                               
## beta_acs_shade[27]      NA                               
## beta_acs_shade[28]      NA                               
## beta_acs_shade[29]      NA                               
## beta_acs_shade[30]      NA                               
## beta_acs_shade[31]      NA                               
## beta_acs_shade[32]      NA                               
## beta_acs_shade[33]      NA                               
## beta_acs_shade[34]      NA                               
## beta_acs_shade[35]      NA                               
## beta_acs_shade[36]      NA                               
## sigma                13.05                               
## sigma_acs             7.57                               
## sigma_shelf           2.98                               
## beta_sp_shade[1]     26.03                               
## beta_sp_shade[2]     14.96                               
## beta_sp_shade[3]     18.69                               
## beta_sp_shade[4]     25.85                               
## beta_sp_shade[5]     22.44                               
## nobs                  1008
plot(coeftab(m.shade.species.sep.acs.int.shelf.no.alpha,m.shade.species.sep.acs.shelf.no.alpha),pars=c(1,44:56))

drop the species shade interaction

m.shade.species.acs.int.shelf <- map2stan(
  alist(
    totleng ~ dnorm(mu,sigma),
    mu <- alpha + 
      alpha_acs[acs_id] +
      alpha_shelf[shelf_id]  +
      beta_sp[species_id] + 
      beta_shade*shade + 
      beta_acs_shade[acs_id]*shade,
    alpha ~ dnorm(53,20),
    alpha_acs[acs_id] ~ dnorm(0,sigma_acs),
    alpha_shelf[shelf_id] ~ dnorm(0,sigma_shelf),
    beta_sp[species_id] ~ dnorm(0,10),
    beta_shade ~ dnorm(0,20),
    beta_acs_shade[acs_id] ~ dnorm(0,10),
    c(sigma,sigma_acs,sigma_shelf) ~ dcauchy(0,1)
  ),
  data = tomato.small,
  iter = 4000, # Rhat a bit high with defaults
  warmup = 1000,
  chains=4,
  cores=4
)
## 
## SAMPLING FOR MODEL 'totleng ~ dnorm(mu, sigma)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 6e-06 seconds (Warm-up)
##                0.000404 seconds (Sampling)
##                0.00041 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 1200 / 12000 ]
[ 2400 / 12000 ]
[ 3600 / 12000 ]
[ 4800 / 12000 ]
[ 6000 / 12000 ]
[ 7200 / 12000 ]
[ 8400 / 12000 ]
[ 9600 / 12000 ]
[ 10800 / 12000 ]
[ 12000 / 12000 ]
save.image("tomato1212.Rdata")
plot(m.shade.species.acs.int.shelf,ask=FALSE)
## Waiting to draw page 2 of 6

## Waiting to draw page 3 of 6

## Waiting to draw page 4 of 6

## Waiting to draw page 5 of 6

## Waiting to draw page 6 of 6

par(mfrow=c(1,1))

precis(m.shade.species.acs.int.shelf,depth = 2)
##                      Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha               41.50   4.95      33.43      49.17  4629    1
## alpha_acs[1]        -2.65   3.32      -7.92       2.70 12000    1
## alpha_acs[2]        -4.97   3.28     -10.14       0.28 12000    1
## alpha_acs[3]        -4.00   3.12      -8.92       0.99 12000    1
## alpha_acs[4]         3.02   3.13      -1.97       7.90 12000    1
## alpha_acs[5]         8.86   3.18       3.52      13.65 12000    1
## alpha_acs[6]         0.18   3.10      -4.69       5.09 12000    1
## alpha_acs[7]         6.43   3.19       1.13      11.28 12000    1
## alpha_acs[8]         0.88   2.93      -3.79       5.46 12000    1
## alpha_acs[9]         2.83   3.01      -1.78       7.80 12000    1
## alpha_acs[10]        0.39   3.49      -5.07       5.98 12000    1
## alpha_acs[11]        3.07   3.02      -1.69       8.04 12000    1
## alpha_acs[12]       -5.60   3.16     -10.64      -0.62 12000    1
## alpha_acs[13]       -0.32   3.44      -5.78       5.17 12000    1
## alpha_acs[14]        1.88   3.20      -3.07       7.13 12000    1
## alpha_acs[15]       -2.30   4.10      -8.56       4.30 12000    1
## alpha_acs[16]        2.69   2.99      -2.24       7.26 12000    1
## alpha_acs[17]        4.89   3.05       0.21       9.92 12000    1
## alpha_acs[18]       -8.27   3.21     -13.04      -2.76 12000    1
## alpha_acs[19]        6.31   3.24       1.33      11.72 12000    1
## alpha_acs[20]       -3.45   2.98      -8.18       1.34 12000    1
## alpha_acs[21]        0.24   2.91      -4.53       4.75 12000    1
## alpha_acs[22]        1.56   4.02      -4.76       7.96 12000    1
## alpha_acs[23]       -1.91   3.51      -7.41       3.63 12000    1
## alpha_acs[24]       -2.62   3.04      -7.48       2.21 12000    1
## alpha_acs[25]       -2.02   2.98      -6.91       2.58 12000    1
## alpha_acs[26]       -1.83   3.35      -7.12       3.59 12000    1
## alpha_acs[27]        1.13   2.98      -3.69       5.84 12000    1
## alpha_acs[28]       -0.90   3.13      -5.82       4.15 12000    1
## alpha_acs[29]       -0.36   3.23      -5.32       4.90 12000    1
## alpha_acs[30]        4.25   3.19      -0.75       9.44 12000    1
## alpha_acs[31]        6.29   3.27       1.19      11.51 12000    1
## alpha_acs[32]       -0.07   3.25      -5.01       5.35 12000    1
## alpha_acs[33]        0.04   3.30      -5.05       5.43 12000    1
## alpha_acs[34]       -7.22   3.18     -12.35      -2.27 12000    1
## alpha_acs[35]       -4.39   3.19      -9.50       0.57 12000    1
## alpha_acs[36]       -2.68   3.97      -8.92       3.69 12000    1
## alpha_shelf[1]      -3.33   2.34      -6.92       0.08  6206    1
## alpha_shelf[2]       4.00   2.38       0.27       7.46  6147    1
## alpha_shelf[3]       0.13   2.34      -3.46       3.57  6271    1
## alpha_shelf[4]      -2.40   2.38      -5.97       1.15  5171    1
## alpha_shelf[5]      -1.34   2.35      -4.91       2.20  5352    1
## alpha_shelf[6]       2.54   2.33      -1.02       6.05  5418    1
## beta_sp[1]          -3.61   4.80     -11.19       4.04  5274    1
## beta_sp[2]           2.37   4.81      -5.23      10.18  5433    1
## beta_sp[3]           3.01   4.77      -4.67      10.50  5439    1
## beta_sp[4]         -11.41   4.87     -19.16      -3.72  5825    1
## beta_sp[5]           6.84   4.80      -0.95      14.49  5381    1
## beta_shade          22.08   3.61      16.73      27.73  4417    1
## beta_acs_shade[1]  -13.40   4.85     -21.04      -5.57 12000    1
## beta_acs_shade[2]  -10.71   4.51     -18.12      -3.71 12000    1
## beta_acs_shade[3]   -3.19   4.06      -9.81       3.06 12000    1
## beta_acs_shade[4]   -5.14   4.14     -11.73       1.48 12000    1
## beta_acs_shade[5]    3.22   4.13      -3.37       9.77 12000    1
## beta_acs_shade[6]    0.45   4.20      -6.38       6.90 12000    1
## beta_acs_shade[7]    6.54   4.17       0.13      13.38 12000    1
## beta_acs_shade[8]   11.25   4.08       5.21      18.24 12000    1
## beta_acs_shade[9]   -1.50   4.07      -8.03       5.03 12000    1
## beta_acs_shade[10]  10.56   5.13       1.99      18.44 12000    1
## beta_acs_shade[11]  -0.61   4.03      -7.27       5.60 12000    1
## beta_acs_shade[12] -12.25   4.53     -19.34      -4.94 12000    1
## beta_acs_shade[13]  -2.57   4.47      -9.48       4.71 12000    1
## beta_acs_shade[14] -10.53   4.73     -17.90      -2.78 12000    1
## beta_acs_shade[15] -11.97   7.18     -23.52      -0.46 12000    1
## beta_acs_shade[16]  -1.49   3.92      -7.62       4.79 12000    1
## beta_acs_shade[17]   3.94   4.22      -2.74      10.64 12000    1
## beta_acs_shade[18]  -2.93   4.44     -10.30       3.83 12000    1
## beta_acs_shade[19]  17.64   4.27      11.04      24.70 12000    1
## beta_acs_shade[20]  -5.75   4.04     -12.21       0.69 12000    1
## beta_acs_shade[21]   0.53   4.03      -5.87       6.98 12000    1
## beta_acs_shade[22]  10.26   5.90       0.82      19.52 12000    1
## beta_acs_shade[23]   4.65   4.76      -2.88      12.43 12000    1
## beta_acs_shade[24]  -7.81   4.03     -14.56      -1.69 12000    1
## beta_acs_shade[25] -10.29   3.93     -16.37      -3.85 12000    1
## beta_acs_shade[26]   8.61   4.55       1.56      15.98 12000    1
## beta_acs_shade[27]   6.76   4.05       0.44      13.24 12000    1
## beta_acs_shade[28]  -3.36   4.25     -10.21       3.37 12000    1
## beta_acs_shade[29]   8.40   4.46       1.30      15.58 12000    1
## beta_acs_shade[30]   8.65   4.10       2.22      15.29 12000    1
## beta_acs_shade[31]   5.16   4.40      -1.94      12.14 12000    1
## beta_acs_shade[32]   8.15   4.54       1.15      15.56 12000    1
## beta_acs_shade[33]  -0.57   4.29      -7.48       6.23 12000    1
## beta_acs_shade[34]  -0.70   4.28      -7.41       6.21 12000    1
## beta_acs_shade[35]  -1.82   4.58      -9.50       5.04 12000    1
## beta_acs_shade[36]  -1.76   5.70     -11.01       7.13 12000    1
## sigma               12.72   0.29      12.28      13.20 12000    1
## sigma_acs            5.06   0.94       3.55       6.46  5823    1
## sigma_shelf          3.54   1.57       1.58       5.43  4851    1
compare(m.shade.species.int.acs.shelf,m.shade.species.acs.int.shelf)
##                                 WAIC pWAIC dWAIC weight    SE   dSE
## m.shade.species.acs.int.shelf 8055.5  63.3   0.0      1 58.42    NA
## m.shade.species.int.acs.shelf 8083.3  42.3  27.8      0 56.39 17.27
coeftab(m.shade.species.acs.int.shelf)
##                    m.shade.species.acs.int.shelf
## alpha                 41.5                      
## alpha_acs[1]         -2.65                      
## alpha_acs[2]         -4.97                      
## alpha_acs[3]            -4                      
## alpha_acs[4]          3.02                      
## alpha_acs[5]          8.86                      
## alpha_acs[6]          0.18                      
## alpha_acs[7]          6.43                      
## alpha_acs[8]          0.88                      
## alpha_acs[9]          2.83                      
## alpha_acs[10]         0.39                      
## alpha_acs[11]         3.07                      
## alpha_acs[12]         -5.6                      
## alpha_acs[13]        -0.32                      
## alpha_acs[14]         1.88                      
## alpha_acs[15]         -2.3                      
## alpha_acs[16]         2.69                      
## alpha_acs[17]         4.89                      
## alpha_acs[18]        -8.27                      
## alpha_acs[19]         6.31                      
## alpha_acs[20]        -3.45                      
## alpha_acs[21]         0.24                      
## alpha_acs[22]         1.56                      
## alpha_acs[23]        -1.91                      
## alpha_acs[24]        -2.62                      
## alpha_acs[25]        -2.02                      
## alpha_acs[26]        -1.83                      
## alpha_acs[27]         1.13                      
## alpha_acs[28]         -0.9                      
## alpha_acs[29]        -0.36                      
## alpha_acs[30]         4.25                      
## alpha_acs[31]         6.29                      
## alpha_acs[32]        -0.07                      
## alpha_acs[33]         0.04                      
## alpha_acs[34]        -7.22                      
## alpha_acs[35]        -4.39                      
## alpha_acs[36]        -2.68                      
## alpha_shelf[1]       -3.33                      
## alpha_shelf[2]           4                      
## alpha_shelf[3]        0.13                      
## alpha_shelf[4]        -2.4                      
## alpha_shelf[5]       -1.34                      
## alpha_shelf[6]        2.54                      
## beta_sp[1]           -3.61                      
## beta_sp[2]            2.37                      
## beta_sp[3]            3.01                      
## beta_sp[4]          -11.41                      
## beta_sp[5]            6.84                      
## beta_shade           22.08                      
## beta_acs_shade[1]    -13.4                      
## beta_acs_shade[2]   -10.71                      
## beta_acs_shade[3]    -3.19                      
## beta_acs_shade[4]    -5.14                      
## beta_acs_shade[5]     3.22                      
## beta_acs_shade[6]     0.45                      
## beta_acs_shade[7]     6.54                      
## beta_acs_shade[8]    11.25                      
## beta_acs_shade[9]     -1.5                      
## beta_acs_shade[10]   10.56                      
## beta_acs_shade[11]   -0.61                      
## beta_acs_shade[12]  -12.25                      
## beta_acs_shade[13]   -2.57                      
## beta_acs_shade[14]  -10.53                      
## beta_acs_shade[15]  -11.97                      
## beta_acs_shade[16]   -1.49                      
## beta_acs_shade[17]    3.94                      
## beta_acs_shade[18]   -2.93                      
## beta_acs_shade[19]   17.64                      
## beta_acs_shade[20]   -5.75                      
## beta_acs_shade[21]    0.53                      
## beta_acs_shade[22]   10.26                      
## beta_acs_shade[23]    4.65                      
## beta_acs_shade[24]   -7.81                      
## beta_acs_shade[25]  -10.29                      
## beta_acs_shade[26]    8.61                      
## beta_acs_shade[27]    6.76                      
## beta_acs_shade[28]   -3.36                      
## beta_acs_shade[29]     8.4                      
## beta_acs_shade[30]    8.65                      
## beta_acs_shade[31]    5.16                      
## beta_acs_shade[32]    8.15                      
## beta_acs_shade[33]   -0.57                      
## beta_acs_shade[34]    -0.7                      
## beta_acs_shade[35]   -1.82                      
## beta_acs_shade[36]   -1.76                      
## sigma                12.72                      
## sigma_acs             5.06                      
## sigma_shelf           3.54                      
## nobs                  1008
plot(coeftab(m.shade.species.acs.int.shelf),pars=c(1:42))

plot(coeftab(m.shade.species.acs.int.shelf),pars=c(43:86))

We don’t need both a species response and an accession response (interaction) term. If we are only keeping one of the terms then the accession interaction model is preferred.

From this I conclude that variation is better considered at the accession level than at the species level.

Smoking deaths among doctors

In 1961 Doll and Hill sent out a questionnaire to all men on the British Medical Register inquiring about their smoking habits. Almost 70% of such men replied. Death certificates were obtained for medical practitioners and causes of death were assigned on the basis of these certificates. The breslow data set contains the person-years of observations and deaths from coronary artery disease accumulated during the first ten years of the study.

Analyse this data set to determine the posterior probability that smoking increases death by coronary artery disease, that age increases death by coronary artery disease, and that there is an interaction between age and smoking.

You can load the data set and learn about the columns using the commands below

data("breslow",package = "boot")
help("breslow",package ="boot")
breslow
##    age smoke     n   y    ns
## 1   40     0 18790   2     0
## 2   50     0 10673  12     0
## 3   60     0  5710  28     0
## 4   70     0  2585  28     0
## 5   80     0  1462  31     0
## 6   40     1 52407  32 52407
## 7   50     1 43248 104 43248
## 8   60     1 28612 206 28612
## 9   70     1 12663 186 12663
## 10  80     1  5317 102  5317

You can think of “person-years” as the number of observations

breslow$n <- as.integer(breslow$n)
m.smoke <- map2stan(alist(
  y~dbinom(n,p),
  logit(p) <- alpha + beta_smoke * smoke,
  alpha ~ dnorm(0,5),
  beta_smoke ~ dnorm(0,5)),
  data=breslow,
  chains = 4,
  cores = 2)
## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used

## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used
## 
## SAMPLING FOR MODEL 'y ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 4e-06 seconds (Warm-up)
##                3.6e-05 seconds (Sampling)
##                4e-05 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
save.image("tomato1212.Rdata")
plot(m.smoke)
pairs(m.smoke)

precis(m.smoke)
##             Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha      -5.96    0.1      -6.12      -5.81  1059    1
## beta_smoke  0.55    0.1       0.39       0.71  1071    1
logistic(-5.97)
## [1] 0.002547734
logistic(-5.97 + 0.55)
## [1] 0.004407633
breslow$age_id <- coerce_index(breslow$age)
m.age <- map2stan(alist(
  y~dbinom(n,p),
  logit(p) <- alpha_age[age_id],
  alpha_age[age_id] ~ dnorm(0,5)),
  data=breslow,
  chains = 4,
  cores = 2)
## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used

## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used
## 
## SAMPLING FOR MODEL 'y ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 4e-06 seconds (Warm-up)
##                4e-05 seconds (Sampling)
##                4.4e-05 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
save.image("tomato1212.Rdata")
plot(m.age)
precis(m.age,depth=2)
##               Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha_age[1] -7.65   0.16      -7.92      -7.40  3234    1
## alpha_age[2] -6.14   0.09      -6.28      -6.00  3800    1
## alpha_age[3] -4.98   0.06      -5.08      -4.88  3655    1
## alpha_age[4] -4.25   0.07      -4.36      -4.15  3571    1
## alpha_age[5] -3.91   0.08      -4.03      -3.77  2974    1
logistic(precis(m.age,depth=2)@output$Mean)
## [1] 0.0004747746 0.0021465847 0.0068064289 0.0140105972 0.0196180444
pairs(m.age)

compare(m.age,m.smoke)
##           WAIC pWAIC dWAIC weight     SE   dSE
## m.age   8614.5   4.5   0.0      1 268.92    NA
## m.smoke 9495.8   1.9 881.3      0 296.55 62.31
m.smoke.age <- map2stan(
  alist(y~dbinom(n,p),
        logit(p) <- alpha_age[age_id] + beta_smoke*smoke,
        alpha_age[age_id] ~ dnorm(0,5),
        beta_smoke ~ dnorm(0,5)),
  data=breslow,
  chains = 4,
  cores = 2)
## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used

## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used
## 
## SAMPLING FOR MODEL 'y ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 5e-06 seconds (Warm-up)
##                4e-05 seconds (Sampling)
##                4.5e-05 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
save.image("tomato1212.Rdata")
plot(m.smoke.age)

pairs(m.smoke.age)

precis(m.smoke.age,depth=2)
##               Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha_age[1] -7.92   0.19      -8.19      -7.60  1732    1
## alpha_age[2] -6.43   0.13      -6.64      -6.24  1221    1
## alpha_age[3] -5.28   0.11      -5.45      -5.10  1046    1
## alpha_age[4] -4.55   0.11      -4.74      -4.38  1058    1
## alpha_age[5] -4.20   0.12      -4.38      -3.99  1163    1
## beta_smoke    0.35   0.11       0.18       0.52   867    1
compare(m.age,m.smoke,m.smoke.age)
##               WAIC pWAIC dWAIC weight     SE   dSE
## m.smoke.age 8604.7   5.6   0.0   0.99 268.47    NA
## m.age       8614.5   4.5   9.8   0.01 268.92  6.60
## m.smoke     9495.8   1.9 891.1   0.00 296.55 61.06
m.smoke.age.int <- map2stan(
  alist(y~dbinom(n,p),
        logit(p) <- alpha_age[age_id] + 
          beta_smoke*smoke + 
          beta_smoke_age[age_id]*smoke,
        alpha_age[age_id] ~ dnorm(0,5),
        beta_smoke ~ dnorm(0,5),
        beta_smoke_age[age_id] ~ dnorm(0,5)),
  data=breslow,
  chains = 4,
  cores = 2)
## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used

## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used
## 
## SAMPLING FOR MODEL 'y ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 6e-06 seconds (Warm-up)
##                4.2e-05 seconds (Sampling)
##                4.8e-05 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
save.image("tomato1212.Rdata")
plot(m.smoke.age.int)
precis(m.smoke.age.int,depth=2)
##                    Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha_age[1]      -9.19   0.70     -10.24      -8.04  2301 1.00
## alpha_age[2]      -6.82   0.29      -7.29      -6.39  2988 1.00
## alpha_age[3]      -5.32   0.19      -5.65      -5.05  3377 1.00
## alpha_age[4]      -4.53   0.19      -4.83      -4.21  2918 1.00
## alpha_age[5]      -3.85   0.18      -4.14      -3.56  3210 1.00
## beta_smoke         0.53   2.04      -2.79       3.68  1004 1.00
## beta_smoke_age[1]  1.25   2.13      -1.99       4.77  1049 1.01
## beta_smoke_age[2]  0.26   2.06      -3.05       3.52  1018 1.00
## beta_smoke_age[3] -0.13   2.05      -3.34       3.14  1005 1.00
## beta_smoke_age[4] -0.21   2.05      -3.57       2.89  1000 1.00
## beta_smoke_age[5] -0.62   2.05      -4.01       2.50  1008 1.00
compare(m.smoke,m.age,m.smoke.age,m.smoke.age.int)
##                   WAIC pWAIC dWAIC weight     SE   dSE
## m.smoke.age.int 8601.6  10.0   0.0   0.83 268.38    NA
## m.smoke.age     8604.7   5.6   3.1   0.17 268.47  7.01
## m.age           8614.5   4.5  12.9   0.00 268.92  9.49
## m.smoke         9495.8   1.9 894.2   0.00 296.55 61.21

pred <- ensemble(m.smoke.age,m.smoke.age.int)
mu <- apply(pred$link,2,mean)
PI <- apply(pred$link,2,PI)
plot.frame <- data.frame(age=breslow$age,smoke=as.factor(breslow$smoke),mean=mu,low=PI[1,],high=PI[2,])
plot.frame
##    age smoke         mean          low         high
## 1   40     0 0.0001677758 0.0000357960 0.0003888649
## 2   50     0 0.0012077086 0.0006502804 0.0017695528
## 3   60     0 0.0049718544 0.0037737892 0.0063698872
## 4   70     0 0.0108106227 0.0077777626 0.0144388881
## 5   80     0 0.0200786706 0.0136071769 0.0268564114
## 6   40     1 0.0005961793 0.0004272274 0.0007811791
## 7   50     1 0.0023859608 0.0020253972 0.0027634900
## 8   60     1 0.0072170599 0.0064127637 0.0080622621
## 9   70     1 0.0146490377 0.0129864984 0.0164390941
## 10  80     1 0.0195327976 0.0165407160 0.0225883393
pl <- ggplot(plot.frame,aes(x=age,y=mean,ymax=high,ymin=low,fill=smoke))
pl <- pl + geom_bar(position = "dodge",stat = "identity")
pl + geom_errorbar(position = position_dodge(width=0.9),width=0.5)

breslow$age_num <- as.numeric(as.character(breslow$age))
m.smoke.age.int.numeric <- map2stan(
  alist(y~dbinom(n,p),
        logit(p) <- alpha +
          beta_age * age_num +
          beta_smoke*smoke + 
          beta_smoke_age*smoke*age_num,
        alpha ~ dnorm(0,10),
        c(beta_smoke,beta_age) ~ dnorm(0,5),
        beta_smoke_age ~ dnorm(0,1)),
  data=breslow,
  chains = 4,
  cores = 2)
## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used

## Warning in FUN(X[[i]], ...): data with name age is not numeric and not used
## 
## SAMPLING FOR MODEL 'y ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 5e-06 seconds (Warm-up)
##                4.2e-05 seconds (Sampling)
##                4.7e-05 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
save.image("tomato1212.Rdata")
plot(m.smoke.age.int.numeric)
precis(m.smoke.age.int.numeric)
##                  Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha          -12.06   0.54     -12.95     -11.22   407 1.02
## beta_smoke       2.05   0.57       1.11       2.94   395 1.02
## beta_age         0.11   0.01       0.09       0.12   412 1.02
## beta_smoke_age  -0.03   0.01      -0.04      -0.01   404 1.02
compare(m.smoke.age.int, m.smoke.age.int.numeric)
##                           WAIC pWAIC dWAIC weight     SE   dSE
## m.smoke.age.int         8601.6  10.0   0.0      1 268.38    NA
## m.smoke.age.int.numeric 8647.8   3.5  46.2      0 269.84 15.02

Cane Sugar

This data comes from an experiment to measure disease resistance in different varieties of sugar cane.

You can get the data and learn about it with:

data("cane",package="boot")
help("cane",package="boot")
head(cane)
##     n  r  x var block
## 1  87 76 19   1     A
## 2 119  8 14   2     A
## 3  94 74  9   3     A
## 4  95 11 12   4     A
## 5 134  0 12   5     A
## 6  92  0  3   6     A
summary(cane)
##        n                r                x              var      block 
##  Min.   : 29.00   Min.   :  0.00   Min.   : 1.00   1      :  4   A:45  
##  1st Qu.: 85.75   1st Qu.:  3.00   1st Qu.: 7.00   10     :  4   B:45  
##  Median :112.50   Median : 11.50   Median :11.50   11     :  4   C:45  
##  Mean   :118.14   Mean   : 20.26   Mean   :11.94   12     :  4   D:45  
##  3rd Qu.:144.50   3rd Qu.: 25.25   3rd Qu.:15.00   13     :  4         
##  Max.   :243.00   Max.   :131.00   Max.   :36.00   14     :  4         
##                                                    (Other):156

Is there evidence of differences in disease resistance in the different varieties? Does including an adaptive prior for Block improve the model fit?

mean(cane$r/cane$n)
## [1] 0.1741397
logit(mean(cane$r/cane$n))
## [1] -1.556568
cane$var_id <- coerce_index(cane$var)
cane$block_id <- coerce_index(cane$block)
cane$n <- as.integer(cane$n)
m.cane.alpha <- map2stan(
  alist(
    r ~ dbinom(n,p),
    logit(p) <- alpha,
    alpha <- dnorm(-1.55,1)
  ),
  data=cane,
  chains = 4,
  cores = 2
)
## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used
## 
## SAMPLING FOR MODEL 'r ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 4e-06 seconds (Warm-up)
##                0.000141 seconds (Sampling)
##                0.000145 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 400 / 4000 ]
[ 800 / 4000 ]
[ 1200 / 4000 ]
[ 1600 / 4000 ]
[ 2000 / 4000 ]
[ 2400 / 4000 ]
[ 2800 / 4000 ]
[ 3200 / 4000 ]
[ 3600 / 4000 ]
[ 4000 / 4000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
plot(m.cane.alpha)
precis(m.cane.alpha)
##        Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha -1.58   0.02       -1.6      -1.54  3562    1
logistic(-1.58)
## [1] 0.1707955

m.var <- map2stan(
  alist(
    r ~ dbinom(n,p),
    logit(p) <- alpha + alpha_var[var_id],
    alpha <- dnorm(-1.55,1),
    alpha_var[var_id] ~ dnorm(0,5)
  ),
  data=cane,
  chains = 4,
  cores = 2,
  iter=10000,
  warmup = 1000
)
## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used
## 
## SAMPLING FOR MODEL 'r ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 4e-06 seconds (Warm-up)
##                0.000165 seconds (Sampling)
##                0.000169 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 3600 / 36000 ]
[ 7200 / 36000 ]
[ 10800 / 36000 ]
[ 14400 / 36000 ]
[ 18000 / 36000 ]
[ 21600 / 36000 ]
[ 25200 / 36000 ]
[ 28800 / 36000 ]
[ 32400 / 36000 ]
[ 36000 / 36000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
save.image("tomato1212.Rdata")
plot(m.var,ask=FALSE)
## Waiting to draw page 2 of 4

## Waiting to draw page 3 of 4

## Waiting to draw page 4 of 4

precis(m.var,depth=2)
##                Mean StdDev lower 0.89 upper 0.89 n_eff Rhat
## alpha         -1.90   0.58      -2.88      -1.01  1072    1
## alpha_var[1]   3.47   0.59       2.56       4.47  1142    1
## alpha_var[2]   0.14   0.59      -0.79       1.12  1126    1
## alpha_var[3]  -2.87   0.82      -4.14      -1.54  2117    1
## alpha_var[4]   0.99   0.59       0.06       1.97  1128    1
## alpha_var[5]   1.17   0.59       0.24       2.13  1106    1
## alpha_var[6]   1.07   0.58       0.17       2.05  1088    1
## alpha_var[7]  -0.71   0.61      -1.66       0.30  1194    1
## alpha_var[8]   0.43   0.58      -0.50       1.39  1102    1
## alpha_var[9]   0.12   0.59      -0.81       1.08  1110    1
## alpha_var[10] -3.26   0.83      -4.55      -1.92  2057    1
## alpha_var[11]  0.68   0.59      -0.24       1.65  1116    1
## alpha_var[12] -0.12   0.59      -1.01       0.89  1129    1
## alpha_var[13]  1.18   0.58       0.27       2.16  1087    1
## alpha_var[14] -0.05   0.59      -0.99       0.91  1119    1
## alpha_var[15]  1.72   0.58       0.85       2.72  1090    1
## alpha_var[16] -2.58   0.77      -3.80      -1.33  1835    1
## alpha_var[17] -0.25   0.59      -1.19       0.72  1135    1
## alpha_var[18] -1.16   0.63      -2.15      -0.13  1268    1
## alpha_var[19]  0.11   0.59      -0.82       1.11  1142    1
## alpha_var[20] -1.93   0.64      -2.94      -0.88  1331    1
## alpha_var[21] -0.66   0.60      -1.60       0.33  1152    1
## alpha_var[22] -1.04   0.65      -2.09      -0.01  1333    1
## alpha_var[23]  3.22   0.59       2.30       4.20  1132    1
## alpha_var[24] -0.03   0.60      -0.97       0.96  1144    1
## alpha_var[25] -4.75   1.20      -6.59      -2.87  3332    1
## alpha_var[26] -1.12   0.60      -2.07      -0.14  1151    1
## alpha_var[27]  0.38   0.59      -0.55       1.37  1139    1
## alpha_var[28] -0.05   0.59      -0.96       0.94  1119    1
## alpha_var[29] -2.05   0.68      -3.14      -0.97  1433    1
## alpha_var[30] -1.57   0.64      -2.61      -0.56  1297    1
## alpha_var[31]  0.96   0.59       0.00       1.90  1129    1
## alpha_var[32]  0.00   0.59      -0.96       0.94  1111    1
## alpha_var[33] -0.37   0.60      -1.33       0.60  1146    1
## alpha_var[34]  0.43   0.59      -0.47       1.43  1116    1
## alpha_var[35]  1.49   0.58       0.55       2.42  1091    1
## alpha_var[36] -1.24   0.65      -2.34      -0.24  1368    1
## alpha_var[37] -0.30   0.59      -1.22       0.69  1132    1
## alpha_var[38] -1.26   0.63      -2.29      -0.27  1250    1
## alpha_var[39] -2.19   0.66      -3.21      -1.10  1371    1
## alpha_var[40]  1.29   0.58       0.40       2.28  1105    1
## alpha_var[41] -1.85   0.64      -2.87      -0.83  1273    1
## alpha_var[42]  0.81   0.58      -0.11       1.78  1103    1
## alpha_var[43] -0.25   0.60      -1.19       0.73  1161    1
## alpha_var[44]  0.94   0.59       0.00       1.90  1113    1
## alpha_var[45]  1.41   0.59       0.49       2.39  1112    1
compare(m.cane.alpha,m.var)
##                 WAIC pWAIC  dWAIC weight     SE   dSE
## m.var        15759.3  45.8    0.0      1 170.57    NA
## m.cane.alpha 19489.1   1.0 3729.9      0 173.20 115.5

m.var.block <- map2stan(
    alist(
      r ~ dbinom(n,p),
      logit(p) <- alpha + alpha_var[var_id] + alpha_block[block_id],
      alpha <- dnorm(-1.55,1),
      alpha_var[var_id] ~ dnorm(0,5),
      alpha_block[block_id] ~ dnorm(0,sigma_block),
      sigma_block ~ dcauchy(0,1)
    ),
    data=cane,
    chains = 4,
    cores = 2,
    iter=10000,
    warmup=4000
  )
## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used
## Warning in FUN(X[[i]], ...): data with name var is not numeric and not used
## Warning in FUN(X[[i]], ...): data with name block is not numeric and not
## used
## 
## SAMPLING FOR MODEL 'r ~ dbinom(n, p)' NOW (CHAIN 1).
## WARNING: No variance estimation is
##          performed for num_warmup < 20
## 
## 
## Chain 1, Iteration: 1 / 1 [100%]  (Sampling)
##  Elapsed Time: 3e-06 seconds (Warm-up)
##                0.000176 seconds (Sampling)
##                0.000179 seconds (Total)
## Computing WAIC
## Constructing posterior predictions
## [ 2400 / 24000 ]
[ 4800 / 24000 ]
[ 7200 / 24000 ]
[ 9600 / 24000 ]
[ 12000 / 24000 ]
[ 14400 / 24000 ]
[ 16800 / 24000 ]
[ 19200 / 24000 ]
[ 21600 / 24000 ]
[ 24000 / 24000 ]
## Aggregated binomial counts detected. Splitting to 0/1 outcome for WAIC calculation.
save.image("tomato1212.Rdata")